 (ns forma.trends.analysis
   (:use [forma.matrix.utils :only (variance-matrix coll-avg)]
-        [clojure.math.numeric-tower :only (sqrt floor abs)]
+        [clojure.math.numeric-tower :only (sqrt floor abs expt)]
         [forma.trends.filter :only (deseasonalize make-reliable hp-filter)])
   (:require [forma.utils :as utils]
             [incanter.core :as i]
             [incanter.stats :as s]))
 
-;; The first few functions are more general, and could be pulled into,
-;; say, the matrix operation namespace.  For now, we will leave them
-;; here, until we can talk more about what goes into that namespace.
-
-(defn time-trend-cofactor
-  "construct a matrix of cofactors; first column is comprised of ones,
-  second column is a range from 1 to [num-months].  The result is a
-  [num-months]x2 incanter matrix."
-  [num-months]
-  (i/trans (i/bind-rows (repeat num-months 1)
-                        (range 1 (inc num-months)))))
+;; TODO (dan): move thest extra matrix functions somewhere else
+
+(defn is-square?
+  "returns true if input matrix is square, and false otherwise"
+  [mat]
+  (let [[row col] (i/dim mat)]
+    (= row col)))
 
 (defn singular?
-  "Check to see if the supplied matrix `X` is singular. note that `X`
-  has to be square, n x n, where n > 1."
-  [X]
-  (<= (i/det X) 0))
-
-(defn ols-coefficient
-  "extract OLS coefficient from a time-series."
-  [ts]
-  (let [ycol (i/trans [ts])
-        X (time-trend-cofactor (count ts))]
-    (i/sel (i/mmult (variance-matrix X)
-                    (i/trans X)
-                    ycol)
-           1 0)))
-
-(defn windowed-apply 
-  "apply a function [f] to a window along a sequence [xs] of length [window]"
-  [f window xs]
-  (pmap f (partition window 1 xs)))
-
-(defn make-monotonic
-  "move through a collection `coll` picking up the min or max values, depending
-  on the value of `comparator` which can be either `min` or `max`.  This is very
-  similar, I think, to the `reduce` function, but with side effects that print
-  into a vector."
-  [comparator coll]
-  (reduce (fn [acc val]
-            (conj acc
-                  (if (empty? acc)
-                    val
-                    (comparator val (last acc)))))
-          []
-          coll))
+  "returns if square input matrix is singular, and false otherwise"
+  [mat]
+  {:pre [(is-square? mat)]}
+  (<= (i/det mat) 0))
+
+(defn idx
+  "return a list of indices starting with 1 equal to the length of
+  input"
+  [coll]
+  (range (count coll)))
+
+(defn average [lst] (/ (reduce + lst) (count lst)))
+
+(defn moving-average [window lst] (map average (partition window 1 lst)))
+
+(defn ols-trend
+  "returns the OLS trend coefficient from the input vector; an
+  intercept is implicitly assumed"
+  [v]
+  (let [time-step (idx v)]
+    (second (:coefs (s/simple-regression v time-step)))))
+
+(defn windowed-map
+  "maps an input function across a sequence of windows onto a vector
+  `v` of length `window-len` and offset by 1"
+  [f window-len v]
+  (pmap f (partition window-len 1 v)))
+
+(defn min-short-trend
+  "returns the minimum value of piecewise linear trends of length
+  `long-block` on a time series `ts`, after smoothing by moving
+  average of window length `short-block`"
+  [long-block short-block ts]
+  (->> (windowed-map ols-trend long-block ts)
+       (moving-average short-block)
+       (apply min)))
+
+(defn long-stats
+  "returns a list with both the value and t-statistic for the OLS
+  trend coefficient for a time series, conditioning on a variable
+  number of cofactors"
+  [ts & cofactors]
+  (let [time-step (idx ts)
+        X (if (empty? cofactors)
+            (i/matrix time-step)
+            (apply i/bind-columns time-step cofactors))]
+    (try
+      (map second (map (s/linear-model ts X) [:coefs :t-tests]))
+      (catch IllegalArgumentException e))))
+
+(defn expt-residuals
+  "returns a list of residuals from the linear regression of `y` on
+  `X`, raised to `power`"
+  [y X power]
+  (map #(expt % power) (:residuals (s/linear-model y X))))
+
+(defn first-order-conditions
+  "returns a matrix with residual weighted cofactors (incl. constant
+  vector) and deviations from mean; corresponds to first-order
+  conditions $f_t$ in the following reference: Hansen, B. (1992)
+  Testing for Parameter Instability in Linear Models, Journal for
+  Policy Modeling, 14(4), pp. 517-533"
+  [y & x]
+  (let [time-step (idx y)
+        X (if (empty? x)
+            (i/matrix time-step)
+            (apply i/bind-columns time-step x))
+        [resid sq-resid] (map (partial expt-residuals y X) [1 2])]
+    (vector (map * resid X)
+            (map * resid (repeat (count y) 1))
+            (map #(- % (average sq-resid)) sq-resid))))
+
+(defn transpose
+  "returns the transposition of a `coll` of vectors"
+  [coll]
+  (apply (partial map vector) coll))
+
+(defn outer-product
+  "returns a flattened vector of the outer product of a vector and its
+  transpose"
+  [coll]
+  (let [mat (i/matrix coll)]
+    (flatten (i/mmult mat (i/trans mat)))))
+
+(defn element-sum
+  "returns a vector of sums of each respective element in a vector of vectors,
+  i.e., a vector with the first element being the sum of the first
+  elements in each sub-vector."
+  [coll]
+  (apply (partial map +) coll))
+
+(defn hansen-mats
+  "returns the matrices of element-wise sums of (1) the first-order
+  conditions, and (2) the cumulative first-order conditions.  This is
+  only an intermediate step in the calculation of the Hansen (1992)
+  test for parameter instability in linear models.
+  TODO: remove redundancy, keep readability"
+  [y & x]
+  (let [foc (apply first-order-conditions y x)]
+    [(element-sum (map outer-product (transpose foc)))
+     (element-sum (map outer-product (transpose (map i/cumulative-sum foc))))]))
+
+(defn hansen-stat
+  "returns the Hansen (1992) test statistic; number of first-order
+  conditions `num-foc`accounts for the demeaned residuals, intercept,
+  and time-step introduced by `first-order-conditions`"
+  [y & x]
+  (let [num-foc (+ (count x) 3) 
+        [foc cumul-foc] (apply hansen-mats y x)]
+    (i/trace
+     (i/mmult
+      (i/solve (i/matrix (map #(* (count y) %) foc) num-foc))
+      (i/matrix cumul-foc num-foc)))))
+
+(defn trend-stats
+  "return a structured map of the trend characteristics; input the
+  spectral time series (e.g., ndvi) and the length (in intervals) of
+  the long and short blocks for the short-term drop measurement, along
+  with the cofactors used in the long-trend calculations and the break
+  magnitude
+
+  NOTE: This is probably not necessary, but it shows what we're trying
+  to do."
+  [y long-block short-block & x]
+  (let [characteristics (create-struct :short :long :long-tstat :break)
+        [long-val long-tstat] (apply long-stats y x)]
+    (struct characteristics
+            (min-short-trend long-block short-block y)
+            long-val
+            long-tstat
+            (apply hansen-stat y x))))
+
+(defn scaled-vector
+  [scalar coll]
+  (map #(* scalar %) coll))
+
+(defn harmonic-series
+  "return two vectors "
+  [freq coll k]
+  (let [pds (count coll)
+        scalar (/ (* 2 (. Math PI) k) freq)]
+    (transpose (map (juxt i/cos i/sin)
+                    (scaled-vector scalar (idx coll))))))
+
+(defn k-harmonic-matrix
+  [freq k coll]
+  (let [degree-vector (vec (map inc (range k)))]
+    (apply i/bind-columns
+           (apply concat
+                  (map (partial harmonic-series freq coll)
+                       degree-vector)))))
+
+(defn harmonic-seasonal-decomposition
+  "returns a deseasonalized time-series; input frequency of
+  data (e.g., 23 for 16-day intervals), number of harmonic terms (use
+  3 to replicate procedure in the Verbesselt (2010)), and the
+  time-series.
+
+  Reference:
+  Verbesselt, J. et al. (2010) Phenological Change Detection while
+  Accounting for Abrupt and Gradual Trends in Satellite Image Time
+  Series, Remote Sensing of Environment, 114(12), 2970–298"
+  [freq k coll]
+  (let [S (:fitted (s/linear-model coll (k-harmonic-matrix freq k coll)))]
+    (map #(+ (average coll) %)
+         (apply (partial map -) [coll S]))))
 
 ;; ## WHOOPBANG :: Collect the short-term drop associated with a
 ;; ## time-series

               product
               (vector product)))))

(defn outer-product
  " calculate the outer product of two vectors

  Example:
  (outer-product [1 2 3] [1 2 3])
  => [1.0000 2.0000 3.0000
      2.0000 4.0000 6.0000
      3.0000 6.0000 9.0000]
"
  [v1 v2]
  (i/matrix
   (for [e v1] (vec (map (partial * e) v2)))))


(defn lengthening-mat

;; deseasonalize the data, which is reflected below.  It would be a
;; good and feasible (easy) bonus to utilize the reliability data.

;; (defn short-trend
;;   "`short-trend` will find the greatest OLS drop over a timeseries `ts`
;;   given sub-timeseries of length `long-block`.  The drops are smoothed
;;   by a moving average window of length `window`."
;;   ([long-block window ts]
;;      (short-trend ts [] long-block window))
;;   ([long-block window ts reli-ts]
;;      (->> (deseasonalize ts)
;;           (windowed-apply ols-coefficient long-block)
;;           (windowed-apply coll-avg window)
;;           (make-monotonic min))))

(defn collect-short-trend
  "preps the short-term drop for the merge into the other data by separating
