module Numeric.Transform.Fourier.Eigensystem where

import qualified Numeric.Transform.Fourier.FFT as FT
import qualified Data.Complex as Complex
import qualified Data.Array as Array
import Data.Array (Array, Ix, )


(^!) :: Num a => a -> Int -> a
(^!) = (^)

{-# specialize gaussianPlain :: Int -> Float  -> Array Int Float #-}
{-# specialize gaussianPlain :: Int -> Double -> Array Int Double #-}

{- |
For small values of @narrow@ (large width)
it is faster to compute the spectrum of the Gaussian
and apply the Fourier transform on it.
-}
gaussianFourier ::
   (Ix a, Integral a, RealFloat b) =>
   a -- ^ array size
   -> b -- ^ reciprocal of width
   -> Array a b -- ^ X[k]
gaussianFourier n narrow =
   fmap ((/ (narrow*sqrt (fromIntegral n))) . Complex.realPart) $
   FT.rfft $ gaussianPlain n (recip narrow)

gaussianPlain ::
   (Ix a, Integral a, RealFloat b) =>
   a -- ^ array size
   -> b -- ^ reciprocal of width, 1 is the width of the eigenvector
   -> Array a b -- ^ X[k]
gaussianPlain n narrow =
   let nr = fromIntegral n
   in  Array.listArray (0,n-1) $
       map (wrappedGaussian 1e-20 (narrow^!2 * nr * pi)) $
       map (/nr) $
       iterate (1+) 0

{- |
For small @c@ the convergence is slow,
but the result will be close to @recip (sqrt c)@.
-}
wrappedGaussian ::
   (RealFloat b) =>
   b -> b -> b -> b
wrappedGaussian tol c x =
   let xpos = iterate (1+) x
       xneg = tail $ iterate (subtract 1) x
       applyTrunc =
          takeWhile ((>=tol) . abs) .
          map (\xi -> exp (- c*xi*xi))
   in  sum (applyTrunc xpos) +
       sum (applyTrunc xneg)


example0 :: ([Double], [Complex.Complex Double])
example0 =
   let c=0.06; n=13 in (map (\k -> (sqrt (pi/c) *) $ sum $ map (\l -> exp (-pi^!2/(n^!2*c)*(k+n*l)^!2)) [-100..100]) [0..n-1], map (\j -> sum $ map (\k -> exp (((-c*k^!2) Complex.:+ 2*pi*j*k/n))) [-100..100]) [0..n-1])

example1 :: (Double, Double)
example1 =
   let c=0.3 in (sum $ map (\k -> exp(-pi*k^!2/c)) [-100..100], sqrt c * sum (map (\k -> exp(-pi*k^!2*c)) [-100..100]))

exampleAlgebraicDouble :: (Double, Double)
exampleAlgebraicDouble =
   (sum $ map (\k -> exp(-pi*k^!2)) [-100..100], sqrt (4 - sqrt 8) * sum (map (\k -> exp(-pi*2*k^!2)) [-100..100]))

{-
I checked polynomials up to degree 9
and there does not seem to be one,
that has the ratio of the numbers as root.

exampleAlgebraicTriple :: (Double, Double)
exampleAlgebraicTriple =
   (sum $ map (\k -> exp(-pi*k^!2)) [-100..100], sum (map (\k -> exp(-pi*3*k^!2)) [-100..100]))
-}

exampleAlgebraicQuadro :: (Double, Double)
exampleAlgebraicQuadro =
   let a = (2-sqrt(2*sqrt 2)) * (2 + sqrt 2)
       _b = 4 / (2 + sqrt(sqrt 8)) :: Double
   in  (sum $ map (\k -> exp(-pi*k^!2)) [-100..100], a * sum (map (\k -> exp(-pi*4*k^!2)) [-100..100]))

exampleAlgebraicHalf :: (Double, Double)
exampleAlgebraicHalf =
   (sum $ map (\k -> exp(-pi*k^!2)) [-100..100], sqrt (2 - sqrt 2) * sum (map (\k -> exp(-pi/2*k^!2)) [-100..100]))

exampleAlgebraicFourth :: (Double, Double)
exampleAlgebraicFourth =
   let a = (2-sqrt(2*sqrt 2)) * (2+sqrt 2)/2
       _b = 2/(2+sqrt(sqrt 8)) :: Double
   in  (sum $ map (\k -> exp(-pi*k^!2)) [-100..100], a * sum (map (\k -> exp(-pi/4*k^!2)) [-100..100]))

exampleAlgebraic2 :: (Double, Double)
exampleAlgebraic2 =
   let n=2 in (sum $ map (\k -> exp(-pi*n*k^!2)) [-100..100], (1+sqrt 2) * sum (map (\k -> exp(-pi*n*(k+1/n)^!2)) [-100..100]))
{-
a = 1+sqrt 2
a - 1 = sqrt 2
(a - 1)^2 = 2
0 = -1 - 2*a + a^2
-}

exampleAlgebraic3 :: (Double, Double)
exampleAlgebraic3 =
   let n=3 in (sum $ map (\k -> exp(-pi*n*k^!2)) [-100..100], (1+sqrt 3) * sum (map (\k -> exp(-pi*n*(k+1/n)^!2)) [-100..100]))
{-
0 = -2 - 2*a + a^2
-}

exampleAlgebraic4 :: [Double]
exampleAlgebraic4 =
   let n=4 in zipWith (*) [1, 1+sqrt(sqrt 2), (1+sqrt(sqrt 2))/(sqrt(sqrt 2)-1), 1+sqrt(sqrt 2)] $ map (\j -> sum (map (\k -> exp(-pi*n*(k+j/n)^!2)) [-100..100])) [0..n-1]
{-
0 = -1 -  4*a + 6*a^2 -  4*a^3 + a^4
0 =  1 - 12*b + 6*b^2 - 12*b^3 + b^4
-}

exampleAlgebraic5 :: [Double]
exampleAlgebraic5 =
   let n=5; a=1.8743053964841243; b=11.833898536015248 in zipWith (*) [1, a, b, a, 1] $ map (\j -> sum (map (\k -> exp(-pi*n*(k+j/n)^!2)) [-100..100])) [0..n-1]
{-
This must be a linear combination of
(goldenRatio, 1, 0, 0, 1)
(goldenRatio, 0, 1, 1, 0)

with goldenRatio = (1+sqrt 5)/2

0=-4*a^0-4*a^1+26*a^2-14*a^3+a^4
0=-4*b^0-4*b^1+26*b^2-14*b^3+b^4
-}



{-
ToDo:
 - efficient computation of 2-norm of the discrete Gaussian
 - investigate eigenfunctions including Gaussian and Hermite polynomials
 - is there an algebraic ratio between exampleAlgebraicTriple and exampleAlgebraicSixtimes
-}