;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; This code models 3-d geometric anisotropy, leading to kriging
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;                       Rotation and stretching....
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
#|
Here's the rotation we'll be using: Mathematica:

{{Cos[z], -Sin[z], 0}, {Sin[z], Cos[z], 0}, {0, 0, 1}}.{{Cos[y], 
     0, -Sin[y]}, {0, 1, 0}, {Sin[y], 0, Cos[y]}}.{{1, 0, 0}, {0, 
        Cos[x], -Sin[x]}, {0, Sin[x], Cos[x]}} // MatrixForm

{
 {Cos[y] Cos[z], -Cos[z] Sin[x] Sin[y] - Cos[x] Sin[z], -Cos[x] Cos[z] Sin[y] + Sin[x] Sin[z]}, 
 {Cos[y] Sin[z],  Cos[x] Cos[z] - Sin[x] Sin[y] Sin[z], -Cos[z] Sin[x] - Cos[x] Sin[y] Sin[z]}, 
 {       Sin[y],                         Cos[y] Sin[x],                         Cos[x] Cos[y]}
 }
|#
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(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 z-axis (counter-clockwise)
     p  - the angle of rotation about the y-axis (counter-clockwise)
     th - the angle of rotation about the x-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       
       (list (* cp ct) (- (+ (* ctsm sp) cmst)) (- smst (* cmct sp))))
      (* 
       s2
       (list (* cp st) (- cmct (* sp smst)) (- (+ ctsm (* cmst sp)))))
      ;; the third direction remains unscaled: it will be scaled "directly" by
      ;; the range value:
      (combine sp (* cp (list sm cm)))
      )
     )
    )
  )
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;                                move-angle
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun move-angle(p)
  "
Move angle simply shifts an angle to the equivalent between -pi and pi.
"
  (- p (* pi (round (/ p pi))))
  )

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Here we use the goliath data of Boyce:
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(def a (read-matrix-data "/home/longa/www/Xlispstat/examples/goliath.txt"))
(def easting (column 2 a))
(def northing (column 3 a))
;; (def aspect (column 4 a))
;; difference highest elevation from lowest:
(def boyce_random_factor (-  3697 3560.784))
(def elevation (* boyce_random_factor (column 5 a)))

(def all
     (apply
      #'bind-columns
      (select 
       (column-lists a)
       (iseq 7 84)
       )
      )
     )
(def svd-all (svd all))

;; (second svd-all)     
(size svd-all)     

(def weights01 (row 1 (transpose (third svd-all))))
(def factor01 (column 1 (first svd-all)))
(def factor02 (column 2 (first svd-all)))
(def factor03 (column 3 (first svd-all)))
(def factor04 (column 4 (first svd-all)))
;; (def factor05 (column 5 (first svd-all)))
;; (def factor06 (column 6 (first svd-all)))
;; (def factor07 (column 7 (first svd-all)))
;; (def factor08 (column 8 (first svd-all)))
;; (def factor09 (column 9 (first svd-all)))

;; (def species 45)
(def species 29)
;; (def species 9)
;; (def species 2)
(def spn (+ species 7))
(def label (strcat "species " (num-to-string species)))
(def species_data (column spn a))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;                  Parameters I'm fiddling with....
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(load-calc "geostats")
;; (def variogram-type #'cosine-variogram)
(def variogram-type #'gaussian-variogram)
;; (def variogram-type #'exponential-variogram)
;; (def variogram-type #'spherical-variogram)
(def iterations 40)
(def lags 20)
(def grid-size 4)
(def yscalemax .1)
(def uni-guess (list .02 400 .03))
(def species_data factor01)
(def label "Factor 01")

(def n (length species_data))
(def positions 
     (transpose
      (list
       easting northing elevation)))
;;			 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;                                nls-function
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun nls-function(vals)
  "
;; For our trial, nls-function computes the value of an anisotropic 
  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
        nugget
"
  (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))
	 (nugget (elt vals 7))
	 )
    ;; Note the use of global variable: variogram-type
    (mapcar
     (lambda(x)
       (apply variogram-type (list x nugget a r))
       )
     ;; Note the use of global variable: species_variogram_positions
     (mapcar
      #'norm
      (column-lists
       (matmult
	(transformation s1 s2 m p th)
	(apply #'bind-columns species_variogram_positions)
	)
       )
      )
     )
    )
  )

(defun univariate-variogram(vals)
  "
;; Univariate-variogram computes the value of an isotropic variogram. 
;; Parameters are
;;	sill
;;	range
;;      nugget
"
  (let* (
	 (sill (elt vals 0))
	 (range (elt vals 1))
	 (nugget (elt vals 2))
	 )
    ;; Note the use of global variable: variogram-type
    (mapcar
     (lambda(x)
       (apply variogram-type (list  x nugget sill range))
       )
     ;; Note the use of global variable: current_bins
     (mapcar #'first current_bins)
     )
    )
  )
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Create the variogram:
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun create-variogram (positions species_data)
  (let (species_variogram)
    (mapcar
     (lambda(m)
       (let* (
	      (point (elt positions m))
	      (species (elt species_data m))
	      (n (length species_data))
	      )
	 (setq 
	  species_variogram
	  (append 
	   species_variogram
	   (mapcar 
	    (lambda(n) 
	      (let* (
		     (point2 (elt positions n))
		     (species2 (elt species_data n))
		     )
		(list (- point point2) (abs (- species species2)))
		)
	      )
	    (iseq (+ m 1) (- n 1))
	    )
	   )
	  )
	 )
       )
     (butlast (iseq n))
     )
    species_variogram
    )
  )
(def species_variogram (create-variogram positions species_data))
;; (def species_variogram (create-variogram positions species_))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Split up the positions and values
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(setq
 species_variogram_positions
 (mapcar #'first species_variogram)
 species_variogram_values
 (mapcar #'second species_variogram)

 ;; here's the three-d case, including elevation:
 species_variogram_distances
 (mapcar #'norm species_variogram_positions)

 max_x (max species_variogram_distances)
 bin-width (/ max_x lags)

 ;; here's the two-d case, pithing elevation:
 species_variogram_distances_2d
 (mapcar (lambda(x) (norm (combine (select x '(0 1)) 0))) species_variogram_positions)

 max_x_2d (max species_variogram_distances_2d)
 bin-width_2d (/ max_x_2d lags)
 )

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Here I bin the sample variogram to create data to fit an isotropic model:
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun create-bins (species_variogram_distances bin-width max_x)
  (mapcar
   (lambda(bin)
     (let* (
	    (diff (- bin bin-width))
	    (keepers
	     (where
	      species_variogram_distances
	      0
	      (lambda(x y) (and (< x bin) (>= x diff))))
	     )
	    )
       (if keepers
	   (list 
	    (mean (select species_variogram_distances keepers))
	    (mean 
	     (select
	      species_variogram_values
	      keepers
	      )
	     )
	    (length keepers)
	    )
	 (list (+ bin (half diff)) 0 0)
	 )
       )
     )
   (rest (rseq 0 max_x lags))
   )
  )
(def bins 
     (create-bins species_variogram_distances bin-width max_x))
(def bins_2d
     (create-bins species_variogram_distances_2d bin-width_2d max_x_2d))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Univariate variogram estimation:
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(def current_bins bins_2d)
(def nreg
     (nreg-model
      #'univariate-variogram
      (mapcar #'second bins_2d)
      uni-guess
      :weights (mapcar #'third bins_2d)
      :parameter-names '(sill range nugget)
      :count-limit iterations
      )
     )
(setq 
 coefs (abs (send nreg :coef-estimates))
 xs (rseq 
     (min (mapcar #'first bins_2d)) 
     (max (mapcar #'first bins_2d)) 
     100)
 nugget (elt coefs 2)
 sill (elt coefs 0)
 range (elt coefs 1)
 )
(xgobips 
 "factor01_2d.ps"
 0
 0
 1
 (mapcar #'first bins_2d)
 (mapcar #'second bins_2d)
 :glyph-type 2
 :glyph-size 5
 :morelines 
 (list 
  (list 1 0 0 
	xs
	(mapcar (lambda(x) (apply variogram-type (list  x nugget sill range))) xs)
	"lndashed"))
 :points (list 
	  (list (repeat 0 lags) (repeat 1 lags) (repeat 0 lags) ;; colors
		(repeat 3 lags) (repeat 3 lags) ;; glyphs
		(mapcar #'first bins_2d)
		(/ (mapcar #'third bins_2d) (sum (mapcar #'third bins_2d)))
		)
	  )
 :title (strcat "Empirical 2d Variogram: " label)
 :xlabel "Distance"
 :ylabel "Same or different"
 :ytics	(rseq 0 yscalemax 6)
 :xmin 0
 :ymin 0
 :ymax yscalemax
 :x-axis t
 :y-axis t
 )

(def current_bins bins)
(def nreg
     (nreg-model
      #'univariate-variogram
      (mapcar #'second bins)
      uni-guess
      :weights (mapcar #'third bins)
      :parameter-names '(sill range nugget)
      :count-limit iterations
      )
     )

(setq 
 coefs (abs (send nreg :coef-estimates))
 xs (rseq 
     (min (mapcar #'first bins)) 
     (max (mapcar #'first bins)) 
     100)
 nugget (elt coefs 2)
 sill (elt coefs 0)
 range (elt coefs 1)
 )

(xgobips 
 "factor01_3d.ps"
 0
 0
 1
 (mapcar #'first bins)
 (mapcar #'second bins)
 :glyph-type 2
 :glyph-size 5
 :morelines 
 (list 
  (list 1 0 0 
	xs
	(mapcar (lambda(x) (apply variogram-type (list  x nugget sill range))) xs)
	"lndashed"))
 :points (list 
	  (list (repeat 0 lags) (repeat 1 lags) (repeat 0 lags) ;; colors
		(repeat 3 lags) (repeat 3 lags) ;; glyphs
		(mapcar #'first bins)
		(/ (mapcar #'third bins) (sum (mapcar #'third bins)))
		)
	  )
 :title (strcat "Empirical 3d Variogram: " label)
 :xlabel "Distance"
 :ylabel "Same or different"
 :ytics	(rseq 0 yscalemax 6)
 :xmin 0
 :ymin 0
 :ymax yscalemax
 :x-axis t
 :y-axis t
 )

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Let's hope that non-linear regression can find a good fit to the parameters:
;; we'll use the univariate parameters as starting values; perhaps that will help.
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(def guess (list sill range 1 1 0 0 0 nugget))
(def nreg
     (nreg-model
      #'nls-function
      species_variogram_values
      guess
      :parameter-names '(sill range scale1 scale2 xangle yangle zangle nugget)
      :count-limit iterations
      )
     )
(setq
 coefs (send nreg :coef-estimates)
 revised 
 (combine 
  (abs (select coefs '(7 0 1 2 3))) (mapcar #'move-angle (select coefs '(4 5 6))))
 )

(mapcar
 (lambda(name value)
   (format t "~%      ~a: ~,4f" name value)
   )
 '(nugget sill range "scale 1" "scale 2" "x-axis rotation" "y-axis rotation" "z-axis rotation")
 revised
 )
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;  Let's take a look at the rotation and scaling:
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(setq
 orth
 (transformation 1 1 
		 (elt revised 5)
		 (elt revised 6)
		 (elt revised 7)
		 )
 trans
 (transformation (elt revised 3)
		 (elt revised 4)
		 (elt revised 5)
		 (elt revised 6)
		 (elt revised 7)
		 )
 )
(mprint orth)
(mprint trans)
;; we know what the svd of trans is!
;; (def svd (svd trans))
;; (mprint (first svd))
;; (print (second svd))
;; (mprint (third svd))

;; in spite of that, I'll call the inverse in my laziness....
(mprint (inverse trans))
(print (mapcar #'norm (column-lists trans)))
(print (mapcar #'norm (column-lists (inverse trans))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;             Now let's transform space, and replot the variogram:
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(setq
 species_variogram_distances_anisotropic
 (let* (
	(a (elt coefs 0))
	(r (elt coefs 1))
	(s1 (elt coefs 2))
	(s2 (elt coefs 3))
	(m (elt coefs 4))
	(p (elt coefs 5))
	(th (elt coefs 6))
	(nugget (elt coefs 7))
	)
   (mapcar 
    #'norm
    (column-lists
     (matmult
      (transformation s1 s2 m p th)
      (apply #'bind-columns species_variogram_positions)
      )
     )
    )
   )
 max_x_anisotropic (max species_variogram_distances_anisotropic)
 bin-width_anisotropic (/ max_x_anisotropic lags)
 bins_anisotropic
 (mapcar
  (lambda(bin)
    (let* (
	   (diff (- bin bin-width_anisotropic))
	   (keepers
	    (where
	     species_variogram_distances_anisotropic
	     0
	     (lambda(x y) (and (< x bin) (>= x diff))))
	    )
	   )
      (if keepers
	  (list 
	   (mean (select species_variogram_distances_anisotropic keepers))
	   (mean 
	    (select
	     species_variogram_values
	     keepers
	     )
	    )
	   (length keepers)
	   )
	(list (+ bin (half diff)) 0 0)
	)
      )
    )
  (rest (rseq 0 max_x_anisotropic lags))
  )
 nugget (elt coefs 7)
 sill (elt coefs 0)
 range (elt coefs 1)
 xs (rseq 
     (min (mapcar #'first bins_anisotropic)) 
     (max (mapcar #'first bins_anisotropic)) 
     100)
 )

(xgobips 
 "factor01_anisotropic_variogram.ps"
 0
 0
 1
 (mapcar #'first bins_anisotropic)
 (mapcar #'second bins_anisotropic)
 :glyph-type 2
 :glyph-size 5
 :morelines 
 (list 
  (list 1 0 0 
	xs
	(mapcar (lambda(x) (apply variogram-type (list x nugget sill range))) xs)
	"lndashed"))
 :points (list 
	  (list (repeat 0 lags) (repeat 1 lags) (repeat 0 lags) ;; colors
		(repeat 3 lags) (repeat 3 lags) ;; glyphs
		(mapcar #'first bins_anisotropic)
		(/ (mapcar #'third bins_anisotropic) (sum (mapcar #'third bins_anisotropic)))
		)
	  )
 :title (strcat "Empirical Variogram following anisotropic shift: " label)
 :xlabel "Transformed Distance"
 :ylabel "Same or different"
 :ytics	(rseq 0 yscalemax 6)
 :xmin 0
 :ymin 0
 :ymax yscalemax
 :x-axis t
 :y-axis t
 )

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;      Now we'll regress against the binned variogram:
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(setq
 current_bins bins_anisotropic
 nreg
 (nreg-model
  #'univariate-variogram
  (mapcar #'second bins_anisotropic)
  (list sill range nugget)
  :weights (mapcar #'third bins_anisotropic)
  :parameter-names '(sill range nugget)
  :count-limit iterations
  )
 coefs (abs (send nreg :coef-estimates))
 xs (rseq 
     (min (mapcar #'first bins_anisotropic)) 
     (max (mapcar #'first bins_anisotropic)) 
     100)
 nugget (elt coefs 2)
 sill (elt coefs 0)
 range (elt coefs 1)
 )
(xgobips 
 "factor01_3d_binned_variogram_model.ps"
 0
 0
 1
 (mapcar #'first bins_anisotropic)
 (mapcar #'second bins_anisotropic)
 :glyph-type 2
 :glyph-size 5
 :morelines 
 (list 
  (list 1 0 0 
	xs
	(mapcar (lambda(x) (apply variogram-type (list  x nugget sill range))) xs)
	"lndashed"))
 :points (list 
	  (list (repeat 0 lags) (repeat 1 lags) (repeat 0 lags) ;; colors
		(repeat 3 lags) (repeat 3 lags) ;; glyphs
		(mapcar #'first bins_anisotropic)
		(/ (mapcar #'third bins_anisotropic) (sum (mapcar #'third bins_anisotropic)))
		)
	  )
 :title (strcat "Empirical 3d isotropic Variogram: " label)
 :xlabel "Transformed Distance"
 :ylabel "Same or different"
 :ytics	(rseq 0 yscalemax 6)
 :xmin 0
 :ymin 0
 :ymax yscalemax
 :x-axis t
 :y-axis t
 )

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;  Now we'll redo those anisotropics based on the improved guess provided by
;;  the binned variogram estimates:
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(def guess (list sill range 
		 (elt revised 3)
		 (elt revised 4)
		 (elt revised 5)
		 (elt revised 6)
		 (elt revised 7)
		 nugget))
(setq
 nreg
 (nreg-model
  #'nls-function
  species_variogram_values
  guess
  :parameter-names '(sill range scale1 scale2 xangle yangle zangle nugget)
  :count-limit iterations
  )
 coefs (send nreg :coef-estimates)
 revised 
 (combine 
  (abs (select coefs '(7 0 1 2 3))) (mapcar #'move-angle (select coefs '(4 5 6))))
 )
(mapcar
 (lambda(name value)
   (format t "~%      ~a: ~,4f" name value)
   )
 '(nugget sill range "scale 1" "scale 2" "x-axis rotation" "y-axis rotation" "z-axis rotation")
 revised
 )
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;  Let's take a look at the rotation and scaling:
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(setq
 orth
 (transformation 1 1 
		 (elt revised 5)
		 (elt revised 6)
		 (elt revised 7)
		 )
 trans
 (transformation (elt revised 3)
		 (elt revised 4)
		 (elt revised 5)
		 (elt revised 6)
		 (elt revised 7)
		 )
 )
(mprint orth)
(mprint trans)
(mprint (inverse trans))
(print (mapcar #'norm (column-lists trans)))
(print (mapcar #'norm (column-lists (inverse trans))))

(setq
 species_variogram_distances_anisotropic
 (let* (
	(a (elt coefs 0))
	(r (elt coefs 1))
	(s1 (elt coefs 2))
	(s2 (elt coefs 3))
	(m (elt coefs 4))
	(p (elt coefs 5))
	(th (elt coefs 6))
	(nugget (elt coefs 7))
	)
   (mapcar 
    #'norm
    (column-lists
     (matmult
      (transformation s1 s2 m p th)
      (apply #'bind-columns species_variogram_positions)
      )
     )
    )
   )
 max_x_anisotropic (max species_variogram_distances_anisotropic)
 bin-width_anisotropic (/ max_x_anisotropic lags)
 bins_anisotropic
 (mapcar
  (lambda(bin)
    (let* (
	   (diff (- bin bin-width_anisotropic))
	   (keepers
	    (where
	     species_variogram_distances_anisotropic
	     0
	     (lambda(x y) (and (< x bin) (>= x diff))))
	    )
	   )
      (if keepers
	  (list 
	   (mean (select species_variogram_distances_anisotropic keepers))
	   (mean 
	    (select
	     species_variogram_values
	     keepers
	     )
	    )
	   (length keepers)
	   )
	(list (+ bin (half diff)) 0 0)
	)
      )
    )
  (rest (rseq 0 max_x_anisotropic lags))
  )
 nugget (elt coefs 7)
 sill (elt coefs 0)
 range (elt coefs 1)
 xs (rseq 
     (min (mapcar #'first bins_anisotropic)) 
     (max (mapcar #'first bins_anisotropic)) 
     100)
 )
(xgobips 
 "factor01_anisotropic_2nd_pass.ps"
 0
 0
 1
 (mapcar #'first bins_anisotropic)
 (mapcar #'second bins_anisotropic)
 :glyph-type 2
 :glyph-size 5
 :morelines 
 (list 
  (list 1 0 0 
	xs
	(mapcar (lambda(x) (apply variogram-type (list  x nugget sill range))) xs)
	"lndashed"))
 :points (list 
	  (list (repeat 0 lags) (repeat 1 lags) (repeat 0 lags) ;; colors
		(repeat 3 lags) (repeat 3 lags) ;; glyphs
		(mapcar #'first bins_anisotropic)
		(/ (mapcar #'third bins_anisotropic) (sum (mapcar #'third bins_anisotropic)))
		)
	  )
 :title (strcat "Empirical Variogram following anisotropic shift: " label)
 :xlabel "Transformed Distance"
 :ylabel "Same or different"
 :ytics	(rseq 0 yscalemax 6)
 :xmin 0
 :ymin 0
 :ymax yscalemax
 :x-axis t
 :y-axis t
 )

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;  Perchance to krig....
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun anisotropic-variogram(x vals)
  (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))
	 (nugget (elt vals 7))
	 (trans (transformation s1 s2 m p th))
	 )
    (mapcar
     (lambda(x) (apply variogram-type (list x nugget a r)))
     (mapcar
      #'norm
      (column-lists
       (matmult
	trans
	(apply #'bind-columns (mapcar (lambda(y) (- x y)) positions))
	)
       )
      )
     )
    )
  )

(def lhs
     (apply
      #'bind-rows
      (mapcar
       (lambda(x) (anisotropic-variogram x coefs))
       positions
       )
      )
     )
(def lhs
     (bind-columns
      (bind-rows
       lhs
       (repeat 1 (second (size lhs)))
       )
      (combine (repeat 1 (first (size lhs))) 0)
      )
     )

(def decomp (lu-decomp lhs))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;            Now we'll get specific:
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(def xgrid (rseq (min easting) (max easting) grid-size))
(def ygrid (rseq (min northing) (max northing) grid-size))
(def zgrid 
     (outer-product
      xgrid
      ygrid
      (lambda(x y)
	(let*  (
		(vec (list x y))
		(weights
		 (mapcar
		  (lambda(v)
		    (let ((diff (- vec v)))
		      (/ (dot diff diff))
		      ))
		  (weave easting northing)
		  )
		 )
		(sum_total (sum weights))
		)
	  (/ 
	   (dot
	    weights
	    elevation
	    )
	   sum_total)
	  )
	)
      )
     )

(def species_estimates (make-array (size zgrid)))
(def
 estimates
 (mapcar
  (lambda(args)
    (let* (
	   (i (first args))
	   (j (second args))
	   (v (select args '(2 3 4)))
	   (weights
	   (lu-solve
	    decomp
	    (combine (anisotropic-variogram v coefs) 1)
	    )
	   )
	  )
      (setf (aref species_estimates i j) 
	    (dot weights (combine species_data 1))
	    )
      )
    )
  (let (stuff)
    (mapcar
     (lambda(i)
       (let ((x (elt xgrid i)))
	 (mapcar
	  (lambda(j)
	    (setq stuff
		  (append stuff 
			  (list
			   (list 
			    i
			    j
			    x
			    (elt ygrid j)
			    (aref zgrid i j)
			    )
			   ))
		  )
	    )
	  (iseq grid-size)
	  )
	 )
       )
     (iseq grid-size)
     )
    stuff
    )
  )
 )
(load-calc "writeR")
(write-r "species.R" species_estimates 'species xgrid ygrid)
(write-vec-r 
 "rawdata.R" 
 (list
  easting northing species_data
  )
  '(raweasting rawnorthing rawdata)
  )