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
^! :: forall a. Num a => a -> Int -> a
(^!) = forall a b. (Num a, Integral b) => a -> b -> 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 :: forall a b. (Ix a, Integral a, RealFloat b) => a -> b -> Array a b
gaussianFourier a
n b
narrow =
   forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((forall a. Fractional a => a -> a -> a
/ (b
narrowforall a. Num a => a -> a -> a
*forall a. Floating a => a -> a
sqrt (forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n))) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Complex a -> a
Complex.realPart) forall a b. (a -> b) -> a -> b
$
   forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a b -> Array a (Complex b)
FT.rfft forall a b. (a -> b) -> a -> b
$ forall a b. (Ix a, Integral a, RealFloat b) => a -> b -> Array a b
gaussianPlain a
n (forall a. Fractional a => a -> a
recip b
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 :: forall a b. (Ix a, Integral a, RealFloat b) => a -> b -> Array a b
gaussianPlain a
n b
narrow =
   let nr :: b
nr = forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n
   in  forall i e. Ix i => (i, i) -> [e] -> Array i e
Array.listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
1) forall a b. (a -> b) -> a -> b
$
       forall a b. (a -> b) -> [a] -> [b]
map (forall b. RealFloat b => b -> b -> b -> b
wrappedGaussian b
1e-20 (b
narrowforall a. Num a => a -> Int -> a
^!Int
2 forall a. Num a => a -> a -> a
* b
nr forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi)) forall a b. (a -> b) -> a -> b
$
       forall a b. (a -> b) -> [a] -> [b]
map (forall a. Fractional a => a -> a -> a
/b
nr) forall a b. (a -> b) -> a -> b
$
       forall a. (a -> a) -> a -> [a]
iterate (b
1forall a. Num a => a -> a -> a
+) b
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 :: forall b. RealFloat b => b -> b -> b -> b
wrappedGaussian b
tol b
c b
x =
   let xpos :: [b]
xpos = forall a. (a -> a) -> a -> [a]
iterate (b
1forall a. Num a => a -> a -> a
+) b
x
       xneg :: [b]
xneg = forall a. [a] -> [a]
tail forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate (forall a. Num a => a -> a -> a
subtract b
1) b
x
       applyTrunc :: [b] -> [b]
applyTrunc =
          forall a. (a -> Bool) -> [a] -> [a]
takeWhile ((forall a. Ord a => a -> a -> Bool
>=b
tol) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => a -> a
abs) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
          forall a b. (a -> b) -> [a] -> [b]
map (\b
xi -> forall a. Floating a => a -> a
exp (- b
cforall a. Num a => a -> a -> a
*b
xiforall a. Num a => a -> a -> a
*b
xi))
   in  forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([b] -> [b]
applyTrunc [b]
xpos) forall a. Num a => a -> a -> a
+
       forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([b] -> [b]
applyTrunc [b]
xneg)


example0 :: ([Double], [Complex.Complex Double])
example0 :: ([Double], [Complex Double])
example0 =
   let c :: Double
c=Double
0.06; n :: Double
n=Double
13 in (forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> (forall a. Floating a => a -> a
sqrt (forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Double
c) forall a. Num a => a -> a -> a
*) forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Double
l -> forall a. Floating a => a -> a
exp (-forall a. Floating a => a
piforall a. Num a => a -> Int -> a
^!Int
2forall a. Fractional a => a -> a -> a
/(Double
nforall a. Num a => a -> Int -> a
^!Int
2forall a. Num a => a -> a -> a
*Double
c)forall a. Num a => a -> a -> a
*(Double
kforall a. Num a => a -> a -> a
+Double
nforall a. Num a => a -> a -> a
*Double
l)forall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100]) [Double
0..Double
nforall a. Num a => a -> a -> a
-Double
1], forall a b. (a -> b) -> [a] -> [b]
map (\Double
j -> forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp (((-Double
cforall a. Num a => a -> a -> a
*Double
kforall a. Num a => a -> Int -> a
^!Int
2) forall a. a -> a -> Complex a
Complex.:+ Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
jforall a. Num a => a -> a -> a
*Double
kforall a. Fractional a => a -> a -> a
/Double
n))) [-Double
100..Double
100]) [Double
0..Double
nforall a. Num a => a -> a -> a
-Double
1])

example1 :: (Double, Double)
example1 :: (Double, Double)
example1 =
   let c :: Double
c=Double
0.3 in (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
kforall a. Num a => a -> Int -> a
^!Int
2forall a. Fractional a => a -> a -> a
/Double
c)) [-Double
100..Double
100], forall a. Floating a => a -> a
sqrt Double
c forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
kforall a. Num a => a -> Int -> a
^!Int
2forall a. Num a => a -> a -> a
*Double
c)) [-Double
100..Double
100]))

exampleAlgebraicDouble :: (Double, Double)
exampleAlgebraicDouble :: (Double, Double)
exampleAlgebraicDouble =
   (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
kforall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100], forall a. Floating a => a -> a
sqrt (Double
4 forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
sqrt Double
8) forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
2forall a. Num a => a -> a -> a
*Double
kforall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
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 :: (Double, Double)
exampleAlgebraicQuadro =
   let a :: Double
a = (Double
2forall a. Num a => a -> a -> a
-forall a. Floating a => a -> a
sqrt(Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a -> a
sqrt Double
2)) forall a. Num a => a -> a -> a
* (Double
2 forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sqrt Double
2)
       _b :: Double
_b = Double
4 forall a. Fractional a => a -> a -> a
/ (Double
2 forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sqrt(forall a. Floating a => a -> a
sqrt Double
8)) :: Double
   in  (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
kforall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100], Double
a forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
4forall a. Num a => a -> a -> a
*Double
kforall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100]))

exampleAlgebraicHalf :: (Double, Double)
exampleAlgebraicHalf :: (Double, Double)
exampleAlgebraicHalf =
   (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
kforall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100], forall a. Floating a => a -> a
sqrt (Double
2 forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
sqrt Double
2) forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Double
2forall a. Num a => a -> a -> a
*Double
kforall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100]))

exampleAlgebraicFourth :: (Double, Double)
exampleAlgebraicFourth :: (Double, Double)
exampleAlgebraicFourth =
   let a :: Double
a = (Double
2forall a. Num a => a -> a -> a
-forall a. Floating a => a -> a
sqrt(Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a -> a
sqrt Double
2)) forall a. Num a => a -> a -> a
* (Double
2forall a. Num a => a -> a -> a
+forall a. Floating a => a -> a
sqrt Double
2)forall a. Fractional a => a -> a -> a
/Double
2
       _b :: Double
_b = Double
2forall a. Fractional a => a -> a -> a
/(Double
2forall a. Num a => a -> a -> a
+forall a. Floating a => a -> a
sqrt(forall a. Floating a => a -> a
sqrt Double
8)) :: Double
   in  (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
kforall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100], Double
a forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Double
4forall a. Num a => a -> a -> a
*Double
kforall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100]))

exampleAlgebraic2 :: (Double, Double)
exampleAlgebraic2 :: (Double, Double)
exampleAlgebraic2 =
   let n :: Double
n=Double
2 in (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
nforall a. Num a => a -> a -> a
*Double
kforall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100], (Double
1forall a. Num a => a -> a -> a
+forall a. Floating a => a -> a
sqrt Double
2) forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
nforall a. Num a => a -> a -> a
*(Double
kforall a. Num a => a -> a -> a
+Double
1forall a. Fractional a => a -> a -> a
/Double
n)forall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100]))
{-
a = 1+sqrt 2
a - 1 = sqrt 2
(a - 1)^2 = 2
0 = -1 - 2*a + a^2
-}

exampleAlgebraic3 :: (Double, Double)
exampleAlgebraic3 :: (Double, Double)
exampleAlgebraic3 =
   let n :: Double
n=Double
3 in (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
nforall a. Num a => a -> a -> a
*Double
kforall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100], (Double
1forall a. Num a => a -> a -> a
+forall a. Floating a => a -> a
sqrt Double
3) forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
nforall a. Num a => a -> a -> a
*(Double
kforall a. Num a => a -> a -> a
+Double
1forall a. Fractional a => a -> a -> a
/Double
n)forall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100]))
{-
0 = -2 - 2*a + a^2
-}

exampleAlgebraic4 :: [Double]
exampleAlgebraic4 :: [Double]
exampleAlgebraic4 =
   let n :: Double
n=Double
4 in forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) [Double
1, Double
1forall a. Num a => a -> a -> a
+forall a. Floating a => a -> a
sqrt(forall a. Floating a => a -> a
sqrt Double
2), (Double
1forall a. Num a => a -> a -> a
+forall a. Floating a => a -> a
sqrt(forall a. Floating a => a -> a
sqrt Double
2))forall a. Fractional a => a -> a -> a
/(forall a. Floating a => a -> a
sqrt(forall a. Floating a => a -> a
sqrt Double
2)forall a. Num a => a -> a -> a
-Double
1), Double
1forall a. Num a => a -> a -> a
+forall a. Floating a => a -> a
sqrt(forall a. Floating a => a -> a
sqrt Double
2)] forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Double
j -> forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
nforall a. Num a => a -> a -> a
*(Double
kforall a. Num a => a -> a -> a
+Double
jforall a. Fractional a => a -> a -> a
/Double
n)forall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100])) [Double
0..Double
nforall a. Num a => a -> a -> a
-Double
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 :: [Double]
exampleAlgebraic5 =
   let n :: Double
n=Double
5; a :: Double
a=Double
1.8743053964841243; b :: Double
b=Double
11.833898536015248 in forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) [Double
1, Double
a, Double
b, Double
a, Double
1] forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Double
j -> forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (\Double
k -> forall a. Floating a => a -> a
exp(-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
nforall a. Num a => a -> a -> a
*(Double
kforall a. Num a => a -> a -> a
+Double
jforall a. Fractional a => a -> a -> a
/Double
n)forall a. Num a => a -> Int -> a
^!Int
2)) [-Double
100..Double
100])) [Double
0..Double
nforall a. Num a => a -> a -> a
-Double
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
-}