#| splsm.lsp 4 Oct 95 ... Updated for MW and pure shared libs. 28 Aug 91 -> 10 Dec 91 |# (require "utils") ;;;; (format t "Loading code...~a~%" (dyn-load "/home/bob/CODERES/splsm.so")) ;;;; (provide "splsm") ;;; ;;; _______________ SPLSM wrapper ______________________ ;;; (defun SPLSM (x y &key (df 6) (order 0) (sorted nil) grid) "Return spline with deg freedom at x's or on input grid. Optional argument for order determines order of differentiation. Need to sort the input X's." ; ... 4 Oct 95 (let* ((ss (apply #'send splsm-proto :new (if sorted (list x y) (let ((ord (order x)) ) (list (select x ord) (select y ord)) )))) (sy (send ss :find-smooth :df df)) (sy (send ss :smooth)) ) (if (or grid (< 0 order)) ; interpolate (if grid (setf sy (send ss :find-deriv grid order)) (setf sy (send ss :find-deriv x order)) ) sy) )) ;;;; #| (def x (rseq 0 (* 2 pi) 100)) (def y (+ (sin x) (* .2 (normal-rand (length x))))) (def p (plot-points x y)) (send p :size 250 200) (def splsm (send splsm-proto :new x y)) (send splsm :collapse) ; fix for dup x's --- none here, so 25 1's (send splsm :find-knots) ; repeats knots at boundary (send splsm :find-smooth :df 4) (def smth4 (send splsm :smooth)) (send p :add-lines x smth4 :color 'red) (send splsm :find-smooth :df 8) (def smth8 (send splsm :smooth)) (send p :add-lines x smth8 :color 'blue) |# (defproto splsm-proto '(x y w ux mapX uy uw knots nKnots df sy lev coef lambda) '(resnum) ; < class slots> *object* ; < parent > ) (defmeth splsm-proto :X-INPUT () (slot-value 'x)) (defmeth splsm-proto :Y-INPUT () (slot-value 'y)) (defmeth splsm-proto :SMOOTH () (if (slot-value 'sy) (if (= (length (slot-value 'sy)) (length (slot-value 'y))) (slot-value 'sy) (select (slot-value 'sy) (where-in (slot-value 'x) (slot-value 'ux))) ))) (defmeth splsm-proto :RESIDUALS () (when (slot-value 'sy) (- (slot-value 'y) (send self :smooth)) )) (defmeth splsm-proto :RESIDUAL-SS () (when (slot-value 'sy) (sum (^ (send self :residuals) 2)) )) (defmeth splsm-proto :ISNEW (x y &optional w) (if (= (length x) (length y)) (let* ((n (length x)) (range (first (- (last x) (first x)))) (digits (+ 6 (floor (* (log 10) (log range))))) (scale (^ 10 digits)) ; round x to 6 decimal digits (x (/ (round (* x scale)) scale)) (ux (coerce (unique x) 'vector)) (last (1- (length ux))) (mapX (/ (- ux (aref ux 0)) (- (aref ux last) (aref ux 0)))) ) (setf (slot-value 'x) (coerce (float x) 'vector)) ; assume sorted (setf (slot-value 'y) (coerce (float y) 'vector)) (setf (slot-value 'w) (coerce (if w (float w) (repeat (float 1) n)) 'vector)) (setf (slot-value 'ux) (float ux)) (setf (slot-value 'mapX) (float mapX)) ) (format t "Dimensions of data are not equal. ~%") )) (defmeth splsm-proto :COLLAPSE () "Fills uy and uw slots with avgs from tied x's." (let* ((n (length (slot-value 'x))) (nu (length (slot-value 'ux))) (indx (where-in (slot-value 'x) (slot-value 'ux) :origin 1)) (scr1 (float (repeat 0.1 nu))) (yb (float (repeat 0.1 nu))) (wb (float (repeat 0.1 nu))) (scr2 (float (repeat 0.1 n))) (result (call-cfun "suff" n nu indx (slot-value 'x) (slot-value 'y) (slot-value 'w) scr1 yb wb scr2) ) ) (setf result (nthcdr 7 result)) (setf (slot-value 'uy) (coerce (first result) 'vector)) (setf (slot-value 'uw) (coerce (second result) 'vector)) )) (defmeth splsm-proto :FIND-KNOTS(&key (allKnots 0)) "Fills knots and nKnots fields." (let* ((nu (length (slot-value 'mapX))) (knots (float (repeat 0.1 (+ 6 nu)))) (result (call-cfun "sknot" ; prob in unix (slot-value 'mapX) nu knots 1 allKnots) )) (setf result (nthcdr 2 result)) (setf (slot-value 'knots) (coerce (first result) 'vector)) (setf (slot-value 'nKnots) (- (first (second result)) 4)) )) (defmeth splsm-proto :FIND-SMOOTH (&key df lambda) "Input deg of freedom or lambda." (unless (slot-value 'uy) (vPrint "Collapsing...~%") (send self :collapse)) (unless (slot-value 'knots) (vPrint "Finding knots...~%") (send self :find-knots)) (if df (if (eql df (slot-value 'df)) ; given df and matches prior df (setf lambda (slot-value 'lambda)) (setf (slot-value 'df) df))) (let* ((nu (length (slot-value 'mapX))) (nk (slot-value 'nKnots)) (coef (repeat 7.7 nk)) (sy (repeat 7.7 nu)) (lev (repeat 7.7 nu)) (icrit (if lambda (list 1 1) (list 3 0))) (lamb (if lambda (float lambda) (float 1 ))) (df (if lambda (float 5) (float df))) (lRng (float (list 0 1.5 0.001))) (scr (repeat 7.7 (* nk (+ 17 nk)))) (dbls (list (float 1) df (float 1) lamb)) (longs (list nu nk 0 4 1 0)) (result (call-cfun "qsbart" dbls longs (slot-value 'mapX) (slot-value 'uy) (slot-value 'uw) (slot-value 'knots) coef sy lev icrit lRng scr) ) (result (coerce result 'vector))) (setf (slot-value 'coef) (aref result 6)) (setf (slot-value 'sy) (aref result 7)) (setf (slot-value 'lev) (aref result 8)) (unless (slot-value 'lambda) (setf (slot-value 'lambda) (first (last (aref result 0))) )) (vPrint "Sum Lev = ~6,3f ; Lambda = ~10,7f ~%" (sum (slot-value 'lev)) (slot-value 'lambda)) )) (defmeth splsm-proto :REPLACE-Y (y) "Replace old y vector if new is same length as prior." (if (= (length y) (length (slot-value 'y))) (progn (setf (slot-value 'y) (coerce (float y) 'vector)) (setf (slot-value 'coef) nil) (setf (slot-value 'sy) nil) (setf (slot-value 'lev) nil) (if (= (length y) (length (slot-value 'uy))) (setf (slot-value 'uy) (slot-value 'y)) (send self :collapse) ) ) (format t "Cannot replace Y; lengths not matched to x.~%") )) (defmeth splsm-proto :FIND-DERIV (x order) (let* ((nux (length (slot-value 'ux))) (min (aref (slot-value 'ux) 0)) (range (- (aref (slot-value 'ux) (1- nux)) min)) (sx (/ (- x min) range)) (sx (select sx (which (<= 0 sx)))) (sx (select sx (which (<= sx 1)))) (len (length sx)) (der (repeat 7.7 len)) (result (call-cfun "bvalues" len (slot-value 'knots) (slot-value 'coef) (slot-value 'nKnots) sx der order)) ) (coerce (first (nthcdr 5 result)) 'vector) ))