; LSTSQ.LSP   (c)1990, Norman Taylor

; This is the complete source code for Norman taylor's Least Squares routine.
; A small portion of this code was included with the article in CADalyst,
; SEPTEMBER 1990. The code is also available directly from the author.

;============= LSTSQ.LSP ===================== REV: 09 JUL 90 =================
;
;  Uses least-squares technique on data set to fit one of curves:
;
;    Y = sum(j: 0, N) { W(j) * P(j, X) },  1 <= N <= 10,  Xmin <= X <= Xmax
;
;    X = sum(j: 0, N) { W(j) * P(j, Y) },  1 <= N <= 10,  Ymin <= Y <= Ymax
;
;  where:  P(j, Z), W(j)  are as defined below.
;
;  Only POINT entities are valid as data points.
;  Equal spacing of the data points along the independent variable axis is
;  not required.
;  User can select order N of fitted curve, dependent variable (X or Y), and
;  number of segments in POLYLINE curve approximation.
;  The curve equation and standard deviation are reported at routine end.
;
;  For the least squares method, the use of orthogonal polynomials is
;  preferable to the coefficient matrix approach, particularly for higher-
;  order equations.  The computational difficulty arises from the value D of
;  the matrix determinant, approximately:
;
;          D = K**n * H(n)
;
;  where:     K = number of data points
;             n = order of coefficient matrix
;          H(n) = Hilbert determinant of order n
;
;  For example:  H(4) = 1.7E-7,  H(6) = 5.4E-18,  H(8) = 2.7E-33
;
;  Even with large numbers of data points, the values D are small enough that
;  roundoff errors in the calculations become more than highly significant.
;
;  There are many sets of orthogonal polynomials, such as the Chebyshev
;  polynomials which require equally-spaced data points.  The generalized
;  orthogonal polynomials, developed by G. E. Forsythe and used in this
;  routine, do not require equally-spaced data points and are defined:
;
;      P(-1, X) = 0,    P(0, X) = 1
;
;      P(1, X) = (X - A(1)) * P(0, X) - B(0) * P(-1, X)
;             ...
;      P(j, X) = (X - A(j)) * P(j-1, X) - B(j-1) * P(j-2, x)
;
;  where:
;
;      S(j) = sum(k: 0, K) { P(j, X(k)) * P(j, X(k)) }
;
;      A(0) = 1,   B(0) = 0
;
;      A(j+1) = sum(k: 0, K) { X(k) * P(j, X(k)) * P(j, X(k)) } / S(j)
;
;      B(j) = sum(k: 0, K) { X(k) * P(j-1, X(k)) * P(j, X(k)) } / S(j-1)
;
;      W(j) = sum(k: 0, K) { Y(k) * P(j, X(k)) } / S(j)
;
;      Y = sum(j: 0, N) { W(j) * P(j, X) }
;
;        = sum(j: 0, N) { C(j) * X**j }
;
;  Ref - "Digital Computation and Numerical Methods", R. W. Southworth and
;        S. W. DeLeeuw, McGraw-Hill 1965, pp.483-5
;
;  Author - Norman F. Taylor, (604)732-4093
;
;============= LIBRARY FUNCTIONS ==============================

;=== "ModeS" saves values of list A sysvars in list ~MLST
(defun ModeS (A / X)
   (setq ~MLST (mapcar '(lambda (X) (list X (getvar X))) A))
)

;=== "ModeN" sets values of list A sysvars
(defun ModeN (A / X)
   (foreach X A (setvar (car X) (cadr X)))
)

;=== "XGetDist" returns dist from opt'l pt or default D
(defun XGetDist (P0 P D)
  (setq P (strcat "\n** " P " <" (rtos D) ">: "))
  (if P0 (cond ((getdist P0 P)) (D)) (cond ((getdist P)) (D)))
)

;=== "XGetInt" returns input integer or default D if input nil
(defun XGetInt (P D)
   (cond ((getint (strcat "\n** " P " <" (itoa D) ">: "))) (D))
)

;=== "XGetKWord" returns input from list OL or default D
(defun XGetKWord (OL P D / K O X)
   (setq K "" O "")
   (foreach X OL
      (setq
         K (strcat K " " X)
         O (if (/= X D) (strcat O "/" X) O)
      )
   )
   (initget (substr K 2))
   (cond
      ((getkword (strcat "\n** " P (substr O 2) "/<" D ">: ")))
      (D)
   )
)

;=== "XGetPoint" returns point or default D if nil response
(defun XGetPoint (P0 P D)
  (setq P (strcat
             "\n** " P " <" (rtos (car D)) "," (rtos (cadr D))
             (if (caddr D) (strcat "," (rtos (caddr D))) "")
             ">: "
          )
  )
  (if P0
     (cond ((getpoint P0 P)) (D))
     (cond ((getpoint P)) (D))
  )
)

;=== "MyError" is standard error handler while command active
(defun MyError (X)
   (if (/= "Function cancelled" X)
      (princ (strcat "\n** Error: " X))
   )
   (command) (command)
   (command "REDRAW" "UNDO" "E")
   (if UNDOIT (command "U"))
   (ModeN ~MLST)
   (setq *Error* OldErr)
   (princ)
)

;============= PROGRAM-SPECIFIC FUNCTIONS =====================

;=== "Poly" returns value of polynomial
(defun Poly (X / Y N)
   (setq N ORD Y (nth N CL))
   (repeat ORD
      (setq N (1- N) Y (+ (* Y X) (nth N CL)))
   )
)

;=== "SumWeights" calculates weights for ortho polynomial sum
(defun SumWeights ()
   (setq
      OPL '()
      OPL (repeat LNG (setq OPL (cons 0.0 OPL)))
      PL (subst 1.0 0.0 OPL)
      BL (list 0.0)
      TMP (float LNG)
      WL (list (/ (apply '+ YL) TMP))
      AL (list 1.0 (/ (apply '+ XL) TMP))
      J 1
   )
   (repeat ORD
      (setq
         A (nth J AL)
         B (nth (1- J) BL)
         OOPL OPL
         OPL PL
         PL (mapcar
               '(lambda (K L M) (- (* (- K A) L) (* B M)))
               XL OPL OOPL
            )
         XPL (mapcar '* XL PL)
         BL (append BL
               (list (/ (apply '+ (mapcar '* XPL OPL)) TMP))
            )
         TMP (apply '+ (mapcar '* PL PL))
         WL (append WL
               (list (/ (apply '+ (mapcar '* YL PL)) TMP))
            )
         AL (append AL
               (list (/ (apply '+ (mapcar '* XPL PL)) TMP))
            )
         J (1+ J)
      )
   )
)

;=== "SumCoefs" collects coefficients for single polynomial
(defun SumCoefs ()
   (setq
      TMP '()
      TMP (repeat ORD (setq TMP (cons 0.0 TMP)))
      OOPL (cons 0.0 TMP)
      OPL (cons 1.0 TMP)
      CL (cons (nth 0 WL) TMP)
      J 1
   )
   (repeat ORD
      (setq
         A (nth J AL)
         B (nth (1- J) BL)
         PL (mapcar
               '(lambda (K L M) (- K (* A L) (* B M)))
               (cons 0.0 OPL) OPL OOPL
            )
         CL (mapcar
               '(lambda (K L) (+ (* (nth J WL) K) L))
               PL CL
            )
         OOPL OPL
         OPL PL
         J (1+ J)
      )
   )
)

;=== "LstSq" derives least-squares fitted curve for data set
(defun LstSq (/ PTS LNG ORD DEP SEG ORG SF XL YL PL OPL OOPL
                XPL AL BL WL CL A B J K L M TMP)
   (princ "\n** Select data points...")
   (setq
      PTS (ssget)                      ;set of data points
      LNG (if PTS (sslength PTS) 0)
      J 0
   )
   (while (< J LNG)                  ;only POINT's in data set
      (setq TMP (ssname PTS J))
      (if (= (cdr (assoc 0 (entget TMP))) "POINT")
         (setq J (1+ J))
         (progn
            (ssdel TMP PTS)
            (setq LNG (1- LNG))
         )
      )
   )
   (setq ORD 0)                        ;order of fitted curve
   (while (or (< ORD 1) (> ORD 10))
      (setq ORD (XGetInt "Order (1 - 10) of curve" 1))
   )
   (if (> LNG ORD)
      (progn                   ;continue if enough data points
         (setq
            TMP (strcat "Order-" (itoa ORD) " in - ")
            DEP (XGetKWord '("X" "Y") TMP "X")
            TMP "Number (> 1) of curve segments"
            SEG (if (> ORD 1) (max (XGetInt TMP 20) 2) 1)
            ORG (XGetPoint nil "Graph origin point" '(0 0))
            SF (list
                  (XGetDist nil "Graph X scale factor" 1.0)
                  (XGetDist nil "Graph Y scale factor" 1.0)
               )
            XL '()
            YL '()
            J 0
         )
         (repeat LNG            ;lists of X- and Y-coords
            (setq
               A (cdr (assoc 10 (entget (ssname PTS J))))
               TMP (mapcar
                      '(lambda (K L M) (/ (- K L) M))
                      (trans A 0 1) ORG SF
                   )
               XL (cons (car TMP) XL)
               YL (cons (cadr TMP) YL)
               J (1+ J)
            )
         )
         (if (= DEP "Y")               ;swap: Y independent var
            (setq TMP XL XL YL YL TMP)
         )
         (SumWeights)
         (SumCoefs)
         (setq B (car XL) A B)     ;independent variable limits
         (foreach TMP XL
            (setq B (max TMP B) A (min TMP A))
         )
         (setq B (/ (- B A) (float SEG)))     ;plot pline curve
         (command "PLINE")
         (setvar "BLIPMODE" 0)
         (setq UNDOIT T)
         (repeat (1+ SEG)
            (setq TMP (list A (Poly A)) A (+ A B))
            (command
               (mapcar
                  '(lambda (K L M) (+ (* K L) M))
                  (if (= DEP "Y") (reverse TMP) TMP) SF ORG
               )
            )
         )
         (command "")
         (princ "\n** Curve:  ")
         (princ (if (= DEP "Y") "X =" "Y ="))
         (setq J -1)
         (repeat (1+ ORD)
            (setq J (1+ J) TMP (nth J CL))
            (princ
               (strcat
                  "\n**         " (if (< TMP 0) "- " "+ ")
                  (rtos (abs TMP)) " * " DEP "**" (itoa J)
               )
            )
         )
         (princ "\n** Standard deviation: ")
         (setq TMP
            (apply '+
               (mapcar
                  '(lambda (J K) (expt (- (Poly J) K) 2))
                  XL YL
               )
            )
         )
         (princ (rtos (sqrt (/ TMP (float (1- LNG))))))
      )
      (princ
         (strcat
            "\n** Minimum of " (itoa (1+ ORD))
            " data points needed"
         )
      )
   )
)

;=== "c:LstSq" is main function
(defun c:LstSq (/ OldErr UNDOIT)
   (setq OldErr *Error* *Error* MyError)
   (ModeS '("CMDECHO" "BLIPMODE" "OSMODE" "UCSFOLLOW"))
   (setvar "CMDECHO" 0)
   (command "UNDO" "GROUP")
   (ModeN '(("OSMODE" 0) ("UCSFOLLOW" 0)))
   (LstSq)                                 ;do least-squares
   (command "REDRAW" "UNDO" "E")
   (ModeN ~MLST)
   (setq *Error* OldErr)
   (princ)
)

(c:LstSq)    ; execute on load

;============= END OF FILE ====================================================