; some utility functions
(defun truncate-inv (x e)
  "Args:  x e.  For an array x, returns 1/x truncating values <= e max(x)"

  (let*  ((perc (* e (max x)))
          (y (make-array (length x) :initial-element 0))
          (z (/ 1 (select x (which (> x perc))))  )
          )
    (dotimes (i (length z))
      (setf (elt y i) (elt z i))
      )
    y)
  )


(defun constrain (u l x)
  "Args:  u l x.  For arrays u, l, x, applies an interval constraint
        to the elements of x:  l <= x <= u."
  (pmin (pmax x l) u)
  )

; this is the basic prototype for our vsp model.
; it has slots for everything I could think of for 
; the simple trunctated svd problem.  Probably the thing
; to do for cross validation is to make a new prototype that
; inherits from this and add slots and methods for the
; cross validation.

(defproto inverse-proto '(jacobian rhs sig eps svd solution 
				   constrained-solution chisq
				   zmin zmax upper-bound lower-bound
				   num-rows num-cols
				   title response) 
  )

; get the slot accessor methods
(load "accessor-methods")

; get the plotting methods
(load  "plot-methods")


(defmeth inverse-proto :describe ()
  (format t (send self :title))
  (format t "~%~d equations and ~d unknowns ~%" 
	     (send self :num-rows)
	     (send self :num-cols)
	     )
  (format t "assumed data standard deviation ~f (ms)~%" (send self :sig))
  (format t "SVD truncation level ~f ~%" (send self :eps))
  )


; svd related methods

; compute the svd and stick it in the svd slot
(defmeth inverse-proto :compute-svd ()
  (send self :svd (sv-decomp (send self :jacobian)))
  )

; this takes the svd  and returns the 
; truncated pseudo-inverse solution.
(defmeth inverse-proto :truncated-svd-solution ()
  (if (send self :svd) ; don't bother if svd slot is empty
      (let* ((tempsvd (send self :svd))
	     (modeleigs (third tempsvd))
	     (singular-values (second tempsvd))
	     (dataeigs (first tempsvd))
	     (temp (matmult modeleigs 
			    (diagonal (truncate-inv
				       singular-values (send self :eps)))
			    (transpose dataeigs) (send self :rhs))))
	(send self :solution temp)
	)
    )
  )