-- ghci Signal.hs -i../linearalgebra/ ../Useful.hs
module Signal where

{-
  This module could benefit from routines in NumericPrelude.Polynomial
  but then we risk usability in Hugs and Lhs2TeX.
-}

import Data.Complex
import Useful(nest,sieve,sliceVert)
import EigensystemNum

import qualified PowerSeries

import System.Cmd(system)
import System.Exit(ExitCode)
import Data.List(intersperse)

import qualified Prelude
import Prelude hiding ((^))

(^) :: Num a => a -> Int -> a
(^) = (Prelude.^)


superpose :: Num a => [a] -> [a] -> [a]
superpose [] ys = ys
superpose xs [] = xs
superpose (x:xs) (y:ys) = x+y : superpose xs ys

{- Superpose two signals where the second one is delayed.
   Cf. MusicSignal, NumericPrelude.LaurentPolynomial -}
delaySuperpose :: Num a => Int -> [a] -> [a] -> [a]
delaySuperpose 0    x     y = superpose x y
delaySuperpose del  []    y = replicate del 0 ++ y
delaySuperpose del (x:xs) y = x : delaySuperpose (del-1) xs y

delaySuperpose1 :: Num a => [a] -> [a] -> [a]
delaySuperpose1 [] [] = []
delaySuperpose1 (x:xs) ys = x : superpose xs ys
delaySuperpose1 _ _ =
   error "delaySuperpose1: convolve can't call it with arguments [] and (_:_)"

convolve :: Num a => [a] -> [a] -> [a]
convolve [] _ = []
convolve (x:xs) ys = delaySuperpose1 (amplify x ys) (convolve xs ys)

amplify :: Num a => a -> [a] -> [a]
amplify x = map (x*)

bsplineMask :: (Num a, Fractional a) => Int -> Int -> [a]
bsplineMask len order = nest order (convolve (replicate len 1)) [1]
--[1/(fromInteger (len^order))]

alternate :: Num a => [a] -> [a]
alternate = zipWith id (cycle [id,negate])

upSample :: Num a => Int -> [a] -> [a]
upSample k = concat . intersperse (replicate (k-1) 0) . map (:[])

downSample :: Int -> [a] -> [a]
downSample = sieve

--zipWith shortens the lists
ziplExt :: (a -> b -> a) -> [a] -> [b] -> [a]
ziplExt _ x [] = x
ziplExt _ [] _ = error "ziplExt: first list must be the longer one"
ziplExt f (x:xs) (y:ys) = f x y : ziplExt f xs ys

foldlCyclic :: (a -> b -> a) -> [a] -> [b] -> [a]
foldlCyclic f acc = foldl (ziplExt f) acc . sliceVert (length acc)
--foldlCyclic f acc x = foldl (ziplExt f) acc (sliceVert (length acc) x)

ziprExt :: (a -> b -> b) -> [a] -> [b] -> [b]
ziprExt _ [] y = y
ziprExt _ _ [] = error "ziprExt: second list must be the longer one"
ziprExt f (x:xs) (y:ys) = f x y : ziprExt f xs ys

foldrCyclic :: (a -> b -> b) -> [b] -> [a] -> [b]
foldrCyclic f acc = foldr (ziprExt f) acc . sliceVert (length acc)



refine :: (Num a) => Int -> [a] -> [a] -> [a]
refine k mask = convolve mask . upSample k


divCeil :: (Integral a) => a -> a -> a
divCeil x y = negate (div (negate x) y)

kadicBand :: Num a => Int -> [a] -> [[a]]
kadicBand shift filt =
   let size    = divCeil (length filt) (shift-1)
       padFilt = replicate (size-1) 0 ++ filt ++ repeat 0
   in  reverse (take size (map (take size) (iterate (drop shift) padFilt)))

adjoint :: RealFloat a => [Complex a] -> [Complex a]
adjoint = map conjugate . reverse

realTransitionFilter :: Num a => [a] -> [a]
realTransitionFilter x = convolve x (reverse x)

transitionFilter :: RealFloat a => [Complex a] -> [Complex a]
transitionFilter x = convolve x (adjoint x)


{-
estimateSpecRadsqr :: Fractional a => [a] -> a
estimateSpecRadsqr x =
  let p1 = sum x
      p2 = sum (map sqr x)
  in  (p1^2-p2)^2/2+p2^2
-}

estimateSpecRadSqr :: Fractional a => [a] -> a
estimateSpecRadSqr x =
  let y  = foldlCyclic (+) [0,0,0] x
      yu = amplify (recip (sum y)) y
      p2 = sum (map (^2) yu)
  in  3/2*(p2-1/3)^2+1/3



{- spectral factorization -}
{- http://sep.stanford.edu/sep/prof/fgdp/c3/paper_html/node5.html#SECTION00130000000000000000 -}
{- Decompose a polynomial into the form convolve (adjoint y) y -}
{- 'x' contains one half of the self-adjoint polynomial -}
spectralFactor :: Floating a => Int -> [a] -> [a]
spectralFactor n x =
   let (b,b2) = spectralFactorRat n x
   in  amplify (sqrt b2) b

spectralFactorRat :: Fractional a => Int -> [a] -> ([a], a)
spectralFactorRat n x =
   let (a,b2) = nest (n-1) (uncurry (levinsonStep x)) ([1], head x)
   in  (take n (PowerSeries.divide [1] (reverse a)), b2)


levinsonStep :: Fractional a => [a] -> [a] -> a -> ([a], a)
levinsonStep r a v =
   let e  = sum (zipWith (*) (tail r) a)
       ev = e/v
       a' = 0 : a
   in (zipWith (-) a' (amplify ev (reverse a')), v - e*ev)


{- move the display routines to SignalIO -}

showTimeValue :: (Num a, Num b) => a -> b -> String
showTimeValue t x = show t ++ " " ++ show x

showSignal :: (Num a, Num b) => a -> a -> [b] -> String
showSignal t dt =
   unlines . zipWith showTimeValue (iterate (dt+) t)

--writeFile "refinablefunction.dat" (showSignal 0 1 (nest 6 (refine 2 [0.4829629, 0.8365163, 0.2241438, -0.1294095]) [1]))


viewSignal :: (Num a) => [a] -> IO ExitCode
viewSignal x =
  let stem = "refinablefunction"
  in  do
        writeFile (stem++".dat") (showSignal (0::Double) 1 x)
        system ("gnuplot "++stem++".gp")
        system ("gs "++stem++".eps")


{-* Test spectral radius estimations for transition matrices. -}

compareSpecRad :: (Fractional a, Ord a) => [a] -> (a,a)
compareSpecRad x =
   (limes (1.0e-12) (specRadApprox (kadicBand 2 x)) ^ 2,
    estimateSpecRadSqr x)

testSpecRad :: (Double, Double)
testSpecRad =
   compareSpecRad (amplify (recip (sqrt 2))
      [0.48296291314469, 0.83651630373747, 0.22414386804186, -0.12940952255092])