{-# LANGUAGE NoMonomorphismRestriction #-}
{-|
Module      : TimeSeries.Scratch
Copyright   : (C) 2013 Parallel Scientific Labs, LLC
License     : GPL-2
Stability   : experimental
Portability : non-portable

Scratches and notes.

-}
module TimeSeries.Scratch where

import Data.List (sort, unfoldr)
import Data.Word
import Text.Printf
import Data.IntMap (IntMap)
import qualified Data.IntMap as IM
import qualified Data.Sequence as Seq

import TimeSeries.PRG64
import TimeSeries.Window (Window)
import TimeSeries.Routing
import qualified TimeSeries.Correlation as C
import qualified TimeSeries.Window as W

-- --------------------------------------------------------------------------
-- * Notes

-- --------------------------------------------------------------------------
-- ** From: /2.2.3 Random Projection/

-- | Size 2 sketch vector of t constructed from random vectors v1 and
-- v2.
--
-- >>> sketch_of_t
-- [0.30000000000000004,-4.58]
--
sketch_of_t :: [Double]
sketch_of_t = map (dotp t) vs where
  v1 = [0.13,  -0.24, 0.47, -0.19]
  v2 = [-0.25, -0.64, 0.17, -0.89]
  vs = [v1,v2]
  t = [1, 2, 3, 4]

dotp :: Num a => [a] -> [a] -> a
dotp a b = sum $ zipWith (*) a b

-- --------------------------------------------------------------------------
-- ** From: /2.3 Matching Pursuit/
--
-- $ Quote:
--
-- \"Given a vector pool containing /n/ time series V = (v_1,v_2,...,v_n), a
-- target vector v_t, tolerated error \\eps, and approximating vector set V_a =
-- \\theta. Define cos\\theta a_1 = ṽ_t * ṽ_1 as the cosine between v_t and a vector v_1 in
-- V. Here /vector v/ = v / ‖v‖.
--
-- 1. Set i = 1;
--
-- 2. Search the pool V and find the vector v_1 whose |cos\theta_1| with
-- respect to v_t is maximal;
--
-- 3. Compute the residue r = v_t - c_1v_1 where c_1 = (‖v_t‖/‖v_1‖) cos\\theta. V_a =
-- V_a ⋃ {v_1}
--
-- 4. If ‖r‖ < \epsilon terminate, return V_a.
--
-- 5. Else set i = i + 1 and v_t = r, go back to 2.\"
--

-- --------------------------------------------------------------------------
-- ** From: /3.4.1 The Sketch Approach/
--
-- $ Quote from the pdf:
--
-- \" ... Quantitatively, given a point x ∈ Rᵐ, we compute its dot product
-- with /d/ random vectors rᵢ ∈ {1, -1}ᵐ. The first random projection of
-- /x/ is given by:
--
-- * y1 = (x*r_1, x*r_2, ..., x*r_d)
--
-- We compute 2b more such random projections y_1,...,y_2b_{+1}. If /w/ is
-- another point in Rᵐ and z_1,...,z_2b_{+1} are its projections using dot
-- products with the same random vectors then the median of
--
-- * ‖y_1-z_1‖, ‖y_2-z_2‖, ..., ‖y_2b_{+1}-z_2b_{+1}‖
--
-- is a good estimate of ‖x-w‖. It lies within a \theta(1/d) factor of ‖x-w‖
-- with probability 1 - (1/2)ᵇ. ... \"
--
-- And,
--
-- \"...Our approach is to use a /\"structured\"/ random vector. The
-- apparently oxymoronic idea is to form each structured random vector r
-- from the concatenation of nb = sw/nb random vectors: r = s_1,...s_{nb},
-- where each si has length bw. Further each si is either u or -u, and u
-- is a random vector in {1, -1}ᵐ. This choice is determined by a random
-- binary nb-vector b: if b = 1, si = u and if bi = 0, si = -u. The
-- structured approach leads to an asymptonic performance of /O(nb)/
-- integer additions and /O(log bw)/ floating point operations per datum
-- and per random vector...\"
--

vx, vw :: [Double]
vx = take 64 $ cycle [19,2,3,4,5,8,7,5,6,1,3,4,2,1,4]
vw = take 64 $ cycle [2,33,1,0,17,8,1,3,1,9,0]

vx', vw' :: [Double]
vx' = [(vxi-mean vx) / variance_direct vx | vxi <- vx]
vw' = [(vwi-mean vw) / variance_direct vw | vwi <- vw]

norm :: Floating a => [a] -> a
norm as = sqrt $ sum [sq a|a<-as]

-- | Direct distance.
dist_direct :: [Double] -> [Double] -> Double
dist_direct x w = norm $ zipWith (-) x w

-- | Sketch distance.
dist_sketch ::
  Int         -- ^ Sketch size /d/.
  -> [Double] -- ^ First vector.
  -> [Double] -- ^ Second vector.
  -> Double
dist_sketch d x w = median norms / sqrt (fromIntegral d) where
  norms = map norm $ zipWith (zipWith (-)) ys zs
  median xs = sort xs !! (length xs `div` 2)
  -- d random vectors with m elements:
  -- ri = {1, -1}ᵐ, i = [1,2,..,2b+1]
  b = 8
  m = length x
  prg64s = chunks m $ concat $ prg64Bits (prg64Init 0x9834219)
  (ys,zs) = foldr f ([],[]) [0..2*b] where
    f n (yacc,zacc) =
      let y = [dotp x r | r <- rs]
          z = [dotp w r | r <- rs]
          rs = take d $ drop (n*d) prg64s
      in  (y:yacc,z:zacc)

chunks :: Int -> [a] -> [[a]]
chunks sz ys
  | null ys   = []
  | otherwise = let (pre,post) = splitAt sz ys in pre : chunks sz post

-- | d / s, where d is direct distance and s is sketch distance.
direct_div_by_sketch_distance :: Integer -> Int -> Double
direct_div_by_sketch_distance seed d =
  let n = 512
      rs = take (n*2) $ prg64Doubles seed
      (x,w) = splitAt n rs
      direct = dist_direct x w
      sketch = dist_sketch d x w
  in  direct / sketch

-- | Print 'direct_div_by_sketch' with different sketch sizes.
print_dist_ratios :: IO ()
print_dist_ratios = do
  let n = 100
      gs = cycle [1..n]
      ds = concatMap (replicate n) [10,20,40,60,80]
      f g d = unwords
        [ "sketch size =", show d ++ ","
        , show (direct_div_by_sketch_distance g d)
        ]
  mapM_ putStrLn $ zipWith f gs ds

-- --------------------------------------------------------------------------
-- ** From: /3.4.2 Partitioning Sketch Vectors/
--
-- $ Quote from the pdf:
--
-- \"... Note that we use correlation and distance more or less
-- interchangebly because one can be computed from the other once the
-- data is normalized. Pearson correlation is related to Euclidean
-- distance as follows:
--
-- > D^2(x̂,ŷ) = 2(1 - corr(x,y))
--
-- Here x̂ and ŷ are obtained from the raw time sereis by computing:\"
--
-- >      x - avg(x)
-- > x̂ = ———————————–
-- >        var(x)

type DistanceFunc = [Double] -> [Double] -> Double
type CorrFunc = [Double] -> [Double] -> Double

-- | From above:
--
-- >                 D^2(x̂,ŷ)
-- > corr(x,y) = 1 - ———————
-- >                    2
--
corr_from_direct_distance :: CorrFunc
corr_from_direct_distance = corr_from_distance dist_direct

-- | Like 'corr_from_direct_distance', but using 'dist_sketch'.
corr_from_sketch_distance ::
  Int
  -- ^ Sketch size.
  -> CorrFunc
corr_from_sketch_distance d = corr_from_distance (dist_sketch d)

-- | Using normalized input vector X and Y.
corr_from_distance :: DistanceFunc -> CorrFunc
corr_from_distance df x y = 1 - ((sq (df x' y'))/(2*len x)) where
  x' = [(xi-mean_x)/std_x|xi <- x]
  y' = [(yi-mean_y)/std_y|yi <- y]
  mean_x = mean x
  mean_y = mean y
  std_x = stddev x
  std_y = stddev y

direct_diff_with_sketch_corr :: Integer -> Int -> (Double, Double)
direct_diff_with_sketch_corr seed d =
  let n = 64
      rs = take (n*2) $ prg64Doubles seed
      (x,y) = splitAt n rs
      direct1 = corr_from_direct_distance x y
      direct2 = correlation_direct_01 x y
      sketch = corr_from_sketch_distance d x y
  in  (abs (direct1 - direct2), abs (direct1 - sketch))

print_corr_diffs :: IO ()
print_corr_diffs = do
  let n = 128
      gs = cycle [1..n]
      ds = concatMap (replicate n) [10,20,40,60,80]
      f g d = unwords
        [ "d=" ++ show d ++ ","
        , let (d1,d2) = direct_diff_with_sketch_corr g d
              p = printf "%0.22f"
          in  "(" ++ p d1 ++ ", " ++ p d2 ++ ")"
        ]
  mapM_ putStrLn $ zipWith f gs ds

-- --------------------------------------------------------------------------
-- ** From: /3.4.3  Algorithmic Framework/
--
-- $ \"... Suppose we are seeking points within some distance d in the
-- original time series space.
--
-- * Partition each sketch vector s of size N into groups of some size g.
--
-- * The ith group of each sketch vector s is placed in the ith grid
-- structure of dimension g.
--
-- * If two sketch vectors s1 and s2 are within distance c ✕ d in more
-- than a fraction f of the groups, then the corresponding windows are
-- candidate highly correlated windows and will be checked exactly.\"
--

-- --------------------------------------------------------------------------
-- ** From: /3.5 The Issues in Implementation/
--
-- $ \"... data normalization and its sketch within a sliding window are as
-- follows (k = sw):
--
-- >       X_k - avg(X_k)
-- > X̂_k = ——————————————
-- >         var(X_k)
-- >
-- > Xsk^k = X̂_k⋅R_k
--
-- ... difficulty lies in that avg(X_k) and var(X_k) change
-- over each basic window .... we will show that the updating is trivial
-- and the sketch needs to be computed only once.\"
--

-- | Update avg and var in a manner described in section 3.5.
incremental_avg_and_var ::
  Integer     -- ^ Random seed
  -> Int      -- ^ /bw/, basic window size
  -> Int      -- ^ /sw/, sliding window size
  -> [Double] -- ^ Input data
  -> [([Double],Double,Double)]
  -- ^ Sliding window contents, avg(Xₖ), and var(Xₖ) (k = sw).
incremental_avg_and_var seed bw sw xs0 = go st0 xs0 where
  -- Number of basic windows.
  nb = sw `div` bw
  -- sw with Floating a => a type.
  sw' = fromIntegral sw

  -- Using whole random vector, should be replaced with incremental
  -- approach. Unused at the moment, sketch is not appearing in this
  -- function yet.
  _r = whole_random_vector seed nb bw

  -- Initial state.
  --
  -- As for a Haskell value, it looks like four elements. 'bws' is
  -- a fixed sized list containing nb elements of previous basic window.
  --
  st0 = (swin0,bws0,sumx0,sumxsq0)

  -- Initial sliding window, filled with zeros.
  swin0 = replicate sw 0
  -- Initial nb basic windows, each basic window is filled with zeros.
  bws0 = replicate nb (replicate bw 0)
  --         nb-1
  -- Initial  ∑ sum(Xₖʲ) (k = sw). From avobe swin0, this value is 0.
  --         j=0
  sumx0 = 0
  --         nb-1
  -- Initial  ∑ sum(Xₖʲ)² (k = sw). From avobe swin0, this value is 0.
  --         j=0
  sumxsq0 = 0

  -- Main recursion, shall be merged in main loop.
  go _ []    = []
  go st xs = (swin',avg',var) : go st' xs' where
    (swin,bws,sumx,sumxsq) = st
    st' = (swin',bws',sumx',sumxsq')
    avg' = (1/sw') * sumx'
    var = ((1/sw') * sumxsq') - sq avg'
    sumx' = sumx - sum bw0 + sum bwnb
    sumxsq' = sumxsq - sum [sq x|x<-bw0] + sum [sq x|x<-bwnb]
    (bwnb,xs') = splitAt bw xs
    bw0 = last bws
    bws' = bwnb : init bws
    swin' = push bwnb swin

-- | Self descriptive function to show the comparison of naive and
-- incremental /avg/ and /var/.
compare_naive_and_incremental_avg_and_var ::
  Integer
  -> Int
  -> Int
  -> [Double]
  -> IO [Double]
compare_naive_and_incremental_avg_and_var seed bw sw xs0 = do
  let f (swin,iavg,ivar) = do
        let navg = mean swin
            nvar = mean [sq x|x<-swin] - (sq navg)
            --
            -- XXX: Can we recover from accumulated rounding errors
            -- after running the system for certain duration? Or we
            -- don't need to care?
            --
            -- > machine_epsilon_mul_10 = 2.22e-16 * 10
            -- > match =
            -- >   abs (iavg-navg) <= machine_epsilon_mul_10 &&
            -- >   abs (ivar-nvar) <= machine_epsilon_mul_10
            --
            -- Instead of comparing with epsilon value, returning
            -- the difference itself to find max in caller function.
            --
            avg_diff = abs $ iavg-navg
            var_diff = abs $ ivar-nvar
        mapM_ putStrLn
          [ "incremental avg : " ++ show iavg
          , "naive avg       : " ++ show navg
          , "incremental var : " ++ show ivar
          , "naive var       : " ++ show nvar
          , ""
          ]
        return $ avg_diff `max` var_diff
  mapM f $ incremental_avg_and_var seed bw sw xs0

-- | Sample comparison of avg and var, as shown in the function name.
comparison_of_avg_and_var :: IO ()
comparison_of_avg_and_var = do
  let seed = 0x2384792
      bw = 32
      sw = 256
      ins = take (sw*512) $ prg64Doubles (succ seed)
  diffs <- compare_naive_and_incremental_avg_and_var seed bw sw ins
  putStrLn $ "Max diff with naive computation: " ++ show (maximum diffs)

-- --------------------------------------------------------------------------
-- ** From: /3.6.3 Performance Tests/
--
-- $ \"... To make the comparison concrete, we should specify our software
-- architecture a bit more. The multi-dimentional search structure we use
-- is in fact a grid structure. The reason we have rejected more
-- sophisticated structures is that we are asking a radius query: which
-- windows (represented as points) are within a certain distance of a
-- given point? A multi-scale structure such as a quadtree or R-tree
-- would not help in this case. Moreover, the grid structure can be
-- stored densely in a hash table so empty cells take up no space. ..\"
--

-- --------------------------------------------------------------------------
-- ** From: /Appendix B: Structured Random Projection for Sliding Window/
--
-- $ \"Given a sliding window length sw and a basic window length bw,
-- we show how to compute the sketches for the sliding windows starting
-- at each timepoint ... we will use a far more efficient approach
-- based on the convolution between \"structured random vector\" and a
-- data vector of length sw.\"
--
-- And from /3.5: The Issues in Implementation/
--
-- $ \"... To normalize the whole sliding window based on the same
-- average and variance, the sketch of the basic window Xbwⁱ,i=0,1,...,nb
-- will be updated as follows.
--
-- >                          bw-1
-- >         (Xbwⁱ⋅R) - avgsw  ∑  Ri
-- >                          i=0
-- > Xskⁱ = ——————————————————————————
-- >                varsw
--
-- where avg will be updated by removing the oldest basic window and
-- adding the new arrival Xnb, ....\"
--

-- | Random vector rbw:
--
-- * rbw = (r₀,r₁,...,rbw-1)
--
-- and /control vector/ b:
--
-- * b = (b₀,b₁,...,bnb-1)
--
unit_random_vector_and_control_vector ::
  Integer    -- ^ Random seed.
  -> Int        -- ^ /nb/, sliding window size / basic window size.
  -> Int        -- ^ /bw/, basic window size.
  -> RandomVector
unit_random_vector_and_control_vector seed nb bw = RandomVector r b where
  r = take bw $ concat $ prg64Bits (prg64Init seed)
  b = take nb $ concat $ prg64Bits (prg64Init (seed+1))

-- | Random vector for sliding window, built from unit random vector
-- and control vector:
--
-- * r = (rbw⋅b₀,rbw⋅b₁,...,rbw⋅bnb-1)
--
whole_random_vector ::
  Integer     -- ^ Random seed.
  -> Int      -- ^ /nb/, sliding window size / basic window size.
  -> Int      -- ^ /bw/, basic window size.
  -> [Double] -- ^ /r/.
whole_random_vector seed nb bw = zipWith (*) b' r' where
  RandomVector r b = unit_random_vector_and_control_vector seed nb bw
  r' = concat $ replicate nb r
  b' = concatMap (replicate bw) b

-- | Random vector used for structured sketches.
data RandomVector = RandomVector
   { -- | Random values treated as /unit/.
     rvUnit :: [Double]
     -- | Random values treated as /controlling vector/.
   , rvControl :: [Double]
   } deriving (Eq, Show)

seq_one_sketch :: Integer -> Int -> Int -> [Double] -> [([Double], Double)]
seq_one_sketch seed bw sw xs = swin_and_sks (swin0,sks0) xs where

   nb = sw `div` bw
   RandomVector r b = unit_random_vector_and_control_vector seed nb bw

   -- Initial contents of sliding window, filled with zeros.
   swin0 = replicate sw 0

   -- Initial nb sketches (inner product between b not computed yet),
   -- filled with zeros.
   sks0 = replicate nb 0

   -- Keeping the state of sks = v₁,v₂,...vnb, as shown in Figure B.6.
   -- Length of sks is fixed to nb. These sketches will be stored per
   -- input stream, per random vector....
   swin_and_sks (swin,sks) ins = case bwin of
     [] -> []
     _  -> (swin', dotp b sks') : swin_and_sks (swin',sks') ins'
     where
       (bwin,ins') = splitAt bw ins
       swin' = push bwin swin
       sks' = dotp r (reverse bwin) : sks

-- | Print comparisions of single sketch.
--
-- Showing sliding window contents, sketch from structured vector,
-- sketch from direct dot product, diff of the two, and ratio of the two.
--
print_comparisons_of_single_sketch :: IO ()
print_comparisons_of_single_sketch = do
  let seed = 0x12387087
      bw = 4
      sw = 16
      r = whole_random_vector seed nb bw
      ins = [1..1024]
      nb = sw `div` bw
  putStrLn "(sliding window,convolved sketch,direct sketch,diff,ratio)"
  sequence_
    [ print (swin,sk,sk',abs(sk-sk'),sk/sk')
    | (swin,sk) <- seq_one_sketch seed bw sw ins
    , let sk' = dotp r swin
    ]

-- | Scratch of:
-- /Figure B.6: Structured convolution procedure every basic window/.
-- Still has \"Curse of Dimensionality\".
--
-- This time, sketch size is controlled with given argument.
--
-- XXX: Normalization is improper.
--
sequence_of_sketches ::
  Integer                  -- ^ Random seed.
  -> Int                   -- ^ Size of sketch.
  -> Int                   -- ^ /bw/, basic window size.
  -> Int                   -- ^ /sw/, sliding window size.
  -> [Double]              -- ^ Input data stream.
  -> [([Double],[Double])] -- ^ Sliding window contents and sketches.
sequence_of_sketches seed d bw sw xs = swin_and_sks (swin0,sks0) xs where

  -- Number of basic windows.
  nb = sw `div` bw

  -- Random unit vectors and control vectors.
  rs_and_bs :: [RandomVector]
  rs_and_bs =
    [ unit_random_vector_and_control_vector (seed+k') nb bw
    | k <- [1..d], let k' = fromIntegral k ]

  -- Random unit vectors and control vectors, separated.
  rss, _bss :: [[Double]]
  (rss,_bss) = foldr f ([],[]) rs_and_bs where
    f (RandomVector rs bs) (racc,bacc) = (rs:racc,bs:bacc)

  -- Initial contents of sliding window, filled with zeros.
  swin0 :: [Double]
  swin0 = replicate sw 0

  -- Initial d sketches of length nb, inner product with control
  -- vector b not computed yet, filled with zeros.
  sks0 :: [[Double]]
  sks0 = replicate d (replicate nb 0)

  -- Keeping the state of sks = v₁,v₂,...vnb, as shown in Figure B.6.
  -- Length of sks is fixed to nb. These sketches will be stored per
  -- input stream.
  swin_and_sks :: ([Double],[[Double]]) -> [Double] -> [([Double],[Double])]
  swin_and_sks (swin,sks) ins = case bwin of
    [] -> []
    _  -> (swin', sketch_of_basicwindow) : swin_and_sks (swin',sks') ins'
    where
      -- Convolved sketch, naive standardization.
      sketch_of_basicwindow = [(xi-mu)/sigma|xi<-xs1]
      xs1 = zipWith f rs_and_bs sks'
      f (RandomVector _rs bs) sk = dotp bs sk
      mu = mean xs1
      sigma = stddev xs1

      -- Convolved sketch ... not properly standardized.
      -- sketch_of_basicwindow = zipWith f rs_and_bs sks'
      -- f (RandomVector rs bs) sk =
      --   (dotp bs sk - (avgsw * sum whole_random_vec)) / varsw where
      --     whole_random_vec = zipWith (*) bs' rs'
      --     bs' = concatMap (replicate bw) bs
      --     rs' = concat $ replicate nb rs
      --     avgsw = mean swin
      --     varsw = stddev swin

      -- Basic window, remaining inputs, and shifted sliding window.
      bwin, ins', swin' :: [Double]
      (bwin, ins') = splitAt bw ins
      swin' = push bwin swin

      -- Pre sketches, inner product with control vector b not
      -- computed yet. These pre sketches get passed to next call of
      -- swin_and_sks.
      sks' :: [[Double]]
      sks' = map (take nb) $ zipWith (:) (dotps rss (reverse bwin)) sks

      dotps :: [[Double]] -> [Double] -> [Double]
      dotps xss ys = [dotp xs' ys | xs' <- xss]

print_convolved_sketches :: Int -> IO ()
print_convolved_sketches d = mapM_ f $ sequence_of_sketches 0 d 4 16 ins where
  ins = take 128 $ prg64Doubles 0xfafafafafa
  f (a,b) = do
    putStrLn "sliding window: "
    print a
    putStrLn $ "average of sketch:"
    print (mean b)
    putStrLn $ "stddev of sketch:"
    print (stddev b)
    putStrLn ""

-- | Push the contents of first list to second list, returning list
-- having same length as second list.
push :: [a] -> [a] -> [a]
push as bs = reverse as ++ take (length bs - length as) bs

-- | Show ratio of correlation value computed with direct function and
-- computed with standardized convolved sketch distance.
--
print_comparisons_incr :: IO ()
print_comparisons_incr = do
  let seed = 0xabcdef789  -- random seed.
      d = 60              -- sketch size.
      sw = 1024           -- sliding window size.
      bw = 128            -- basic window size.
      nb = sw `div` bw

      -- whole random vectors with using same seeds as in
      -- 'sequence_of_sketches'.
      _rs = [whole_random_vector (seed+k) nb bw | k <- [1..d]]

      -- Input values, to represent single time series.
      ins1 = take 8192 $ prg64Doubles 0x7777777
      epsilons = prg64Doubles 0x8888888
      ins2 = take 8192 $ [x+(e*0.5)|(x,e) <- zip ins1 epsilons]

  sequence_
    [ putStrLn (printf "(direct,sketch): (%.22f, %.22f)" direct sketch) >>
      putStrLn (printf "ratio: %.22f" ratio) >>
      putStrLn (printf "(avg1,var1): (%.22f, %.22f)" avg1 var1) >>
      putStrLn (printf "(avg2,var2): (%.22f, %.22f)" avg2 var2) >>
      putStrLn ""
    | ((swin1s,sk1s),(swin2s,sk2s)) <-
         -- Same seed value is required
         zip (sequence_of_sketches seed d bw sw ins1)
             (sequence_of_sketches seed d bw sw ins2)
    , let direct = C.direct (W.fromList swin1s) (W.fromList swin2s)
    , let sketch = 1 - (d'/len') where
            d' = sum [sq (x-y)|(x,y) <- zip sk1s sk2s]
            len' = 2 * len sk1s
    , let ratio = (direct/sketch) :: Double
    , let (avg1,var1) = (mean sk1s, stddev sk1s)
    , let (avg2,var2) = (mean sk2s, stddev sk2s)
    ]

-- --------------------------------------------------------------------------
-- ** From: /Appendix C: An Upper Bound of the Grid Size/
--
-- $ Quote:
--
-- \"Now let's examine how to embed this sketch filter in a grid
-- structure.  At first, we assume our parameter group is (N,g,c,f)
-- and therfore there are totally ng = N/g groups. We will assign one
-- grid structure to each sketch group. For each grid structure the
-- critical parameter is the largest value A and diameter a. Now let's
-- bound the size of A...\".
--

-- --------------------------------------------------------------------------
-- * Scratches

w0, w1 :: Window
w0 = W.fromList $ take 128 (cycle [1,2,3,4])
w1 = W.empty 32

-- For copying contents of one window to another relatively smaller
-- window.
w1_copied_from_w0 :: Window
w1_copied_from_w0 = W.copyContents w1 w0

len :: Num b => [a] -> b
len zs = fromIntegral (length zs)

mean :: Fractional a => [a] -> a
mean zs = sum zs / len zs

stddev :: [Double] -> Double
stddev zs = sqrt (mean [sq zi|zi<-zs] - (sq (mean zs)))

variance_direct :: [Double] -> Double
variance_direct x = inner / (len x) where
  x_minus_mean = [xi-mean x| xi <- x]
  inner = dotp x_minus_mean x_minus_mean

covariance_direct :: [Double] -> [Double] -> Double
covariance_direct x y = innerp / (len x) where
  innerp = dotp [xi-mean_x|xi<-x] [yi-mean_y|yi<-y]
  mean_x = mean x
  mean_y = mean y

correlation_direct_01 :: [Double] -> [Double] -> Double
correlation_direct_01 x y = cov / sqrt (var_x * var_y) where
  cov = covariance_direct x y
  var_x = variance_direct x
  var_y = variance_direct y

correlation_direct_02 :: [Double] -> [Double] -> Double
correlation_direct_02 x y = dotp (f x) (f y) where
  f z =
    [ z'i / z'_sqrt
    | let mean_z = mean z
    , let z' = [zi-mean_z| zi <- z]
    , let z'_sqrt = sqrt (dotp z' z')
    , z'i <- z' ]

standardize :: [Double] -> [Double]
standardize xs = [xi - mean xs / variance_direct xs | xi <- xs]

-- | Generate list of doubles between 0 to 1 with 'PRG64'.
prg64Doubles :: Integer -> [Double]
prg64Doubles = unfoldr f . prg64Init where
  f seed =
    let (seed',w64) = prg64Next seed
    in  Just (fromIntegral w64 / max64, seed')
  max64 = fromIntegral (maxBound :: Word64)

bwins01 :: BasicWindows
bwins01 = IM.fromList $ zip [0..] wins where
  wins = map (f base) [0,1,2,3]
  f xs k = Seq.fromList [W.fromList [ai+k|ai<-as] | as <- xs]
  base = [[0..9],[10..19],[20..29],[30..39]]

swins01 :: IntMap Window
swins01 = slidingWindow 40 bwins01

sq :: Floating a => a -> a
sq n = n ^ (2::Int)