;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;                       Rotation and stretching....
#|
Not trusting myself to get the result right, I did this in Mathematica:

{{Cos[t], -Sin[t], 0}, {Sin[t], Cos[t], 0}, {0, 0, 1}}.{{Cos[p], 
        0, -Sin[p]}, {0, 1, 0}, {Sin[p], 0, Cos[p]}}.{{1, 0, 0}, {0, 
        Cos[m], -Sin[m]}, {0, Sin[m], Cos[m]}} // MatrixForm

{{Cos[p] Cos[t], -Cos[t] Sin[m] Sin[p] - 
      Cos[m] Sin[t], -Cos[m] Cos[t] Sin[p] + Sin[m] Sin[t]}, {Cos[p] Sin[t], 
    Cos[m] Cos[t] - Sin[m] Sin[p] Sin[t], -Cos[t] Sin[m] - 
      Cos[m] Sin[p] Sin[t]}, {Sin[p], Cos[p] Sin[m], Cos[m] Cos[p]}}
|#
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun transformation(s1 s2 m p th )
  "
;; Transformation computes the transformation matrix to turn an ellipsoid into
      a sphere, so that we can use an isotropic (omni-directional) variogram
      model.

     s1 - scale factor multiplying the x-direction; 1/s1 multiplies y
     s2 - scale factor multiplying the y-direction; 1/s2 multiplies z
     m  - the angle of rotation about the x-axis (counter-clockwise)
     p  - the angle of rotation about the y-axis (counter-clockwise)
     th - the angle of rotation about the z-axis (counter-clockwise)

;; Examples:
 (transformation 1 1 0 0 0)
 (transformation 1 1 pi/4 0 0)
 (transformation 1 1 0 pi/4 0)
 (transformation 1 1 0 0 pi/4)

 (matmult (transformation 1 1 pi/4 0 0) (list 1 1 1))
 (matmult (transformation 1 1 0 pi/4 0) (list 1 1 1))
 (matmult (transformation 1 1 0 0 pi/4) (list 1 1 1))
"
  (let* (
	 (cp (cos p))
	 (ct (cos th))
	 (cm (cos m))
	 (sp (sin p))
	 (st (sin th))
	 (sm (sin m))
	 (cmct (* cm ct))
	 (smst (* sm st))
	 (ctsm (* ct sm))
	 (cmst (* cm st))
	 )
    (matrix
     '(3 3)
     (combine
      (* s1        (combine (* cp ct) (- (+ (* ctsm sp) cmst)) (- smst (* cmct sp))))
      (* (/ s2 s1) (combine (* cp st) (- cmct (* sp smst)) (- (+ ctsm (* cmst sp)))))
      (* (/ s2)    (combine sp (* cp (list sm cm))))
      )
     )
    )
  )

(defun nls-function(vals)
  "
;; For our trial, nls-function computes the value of an isotropic gaussian
  variogram, after transforming space. Parameters are
	a - sill
	r - range
	s1 - scale 1
	s2 - scale 2
	m - x-axis rotation
	p - y-axis rotation
	th - z-axis rotation

;; Examples:
 (def targets (nls-function (list 3 10 1 2 .1 .2 .3)))
"
  (let* (
	 (a (elt vals 0))
	 (r (elt vals 1))
	 (s1 (elt vals 2))
	 (s2 (elt vals 3))
	 (m (elt vals 4))
	 (p (elt vals 5))
	 (th (elt vals 6))
	 (transformation (transformation s1 s2 m p th))
	 )
    (mapcar
     (lambda(x)
       (* a (exp (- (square (/ (norm (matmult transformation x)) r)))))
       )
     positions
     )
    )
  )

;; Here are 100 positions in 3-space, generated using a random normal:
(def positions (mapcar (lambda(x) (rand.n 3)) (iseq 100)))

;; These are the actual parameters used to generate the target:
(def parameters (list 3 10 3 2 .1 .2 .3))
(def targets (nls-function parameters))
;; (so nls-function has transformed the normal cloud of locations according to
;; the parameters, then computes the value of the isotropic gaussian variogram
;; for a target to our non-linear regression procedure.

;; We provide a guess, close but away from the actual parameters above.
;; I'm imagining that the 2 and 8 are coming out of a omni-directional variogram
;; model, and that we're using the "null" starting position of omni-directional
;; for our choices of the scale factors and the angles - so 1 1 0 0 0
(def guess (list 2 8 1 1 0 0 0))

;; and hope that non-linear regression can find the actual parameters: it does!
(nreg-model
 #'nls-function
 targets
 guess
 :parameter-names '(a r s1 s2 xr yr zr)
 )