packages feed

numeric-prelude 0.1 → 0.1.1

raw patch · 18 files changed

+1156/−140 lines, 18 filesdep +HTamdep +gnuplotdep +storable-recordnew-component:exe:test-gaussian

Dependencies added: HTam, gnuplot, storable-record

Files

numeric-prelude.cabal view
@@ -1,5 +1,5 @@ Name:           numeric-prelude-Version:        0.1+Version:        0.1.1 License:        GPL License-File:   LICENSE Author:         Dylan Thurston <dpt@math.harvard.edu>, Henning Thielemann <numericprelude@henning-thielemann.de>, Mikael Johansson@@ -131,9 +131,12 @@   default:     False  Library-  Build-Depends: parsec >=1 && <3, HUnit >=1 && <2, QuickCheck >=1 && <2-  Build-Depends: non-negative >=0.0.2 && <0.1-  Build-Depends: utility-ht >=0.0.4 && <0.1+  Build-Depends:+    parsec >=1 && <3,+    QuickCheck >=1 && <2,+    storable-record >=0.0.1 && <0.1,+    non-negative >=0.0.2 && <0.1,+    utility-ht >=0.0.4 && <0.1   If flag(splitBase)     Build-Depends:       base >= 2 && <5,@@ -222,10 +225,12 @@     Number.Physical     Number.Physical.Read     Number.Physical.Show+    NumericPrelude.Elementwise     NumericPrelude     PreludeBase   Other-modules:     NumericPrelude.List+    Algebra.AffineSpace     MathObj.Gaussian.Variance     MathObj.Gaussian.Bell     MathObj.Gaussian.Polynomial@@ -247,5 +252,19 @@     Test.MathObj.Gaussian.Variance     Test.MathObj.Gaussian.Bell   Main-Is: Test/Run.hs-  If !flag(buildTests)-    Buildable:         False+  If flag(buildTests)+    Build-Depends: HUnit >=1 && <2+  Else+    Buildable: False++Executable test-gaussian+  Hs-Source-Dirs: src, test+  Main-Is: Gaussian.hs+  Other-Modules:+    MathObj.Gaussian.Example+  If flag(buildTests)+    Build-Depends:+      gnuplot >=0.3 && <0.4,+      HTam >=0.0.2 && <0.1+  Else+    Buildable: False
src/Algebra/Additive.hs view
@@ -1,15 +1,19 @@ {-# LANGUAGE NoImplicitPrelude #-} module Algebra.Additive (-    {- * Class -}+    -- * Class     C,     zero,     (+), (-),     negate, subtract, -    {- * Complex functions -}+    -- * Complex functions     sum, sum1, -    {- * Instances for atomic types -}+    -- * Instance definition helpers+    elementAdd, elementSub, elementNeg,+    (<*>.+), (<*>.-), (<*>.-$),++    -- * Instances for atomic types     propAssociative,     propCommutative,     propIdentity,@@ -21,6 +25,10 @@ import Data.Int  (Int,  Int8,  Int16,  Int32,  Int64,  ) import Data.Word (Word, Word8, Word16, Word32, Word64, ) +import qualified NumericPrelude.Elementwise as Elem+import Control.Applicative (Applicative(pure, (<*>)), )+import Data.Tuple.HT (fst3, snd3, thd3, )+ import qualified Data.Ratio as Ratio98 import qualified Prelude as P import Prelude(Int, Integer, Float, Double, fromInteger, )@@ -86,9 +94,63 @@   +{- |+Instead of baking the add operation into the element function,+we could use higher rank types+and pass a generic @uncurry (+)@ to the run function.+We do not do so in order to stay Haskell 98+at least for parts of NumericPrelude.+-}+{-# INLINE elementAdd #-}+elementAdd ::+   (C x) =>+   (v -> x) -> Elem.T (v,v) x+elementAdd f =+   Elem.element (\(x,y) -> f x + f y) -{-* Instances for atomic types -}+{-# INLINE elementSub #-}+elementSub ::+   (C x) =>+   (v -> x) -> Elem.T (v,v) x+elementSub f =+   Elem.element (\(x,y) -> f x - f y) +{-# INLINE elementNeg #-}+elementNeg ::+   (C x) =>+   (v -> x) -> Elem.T v x+elementNeg f =+   Elem.element (negate . f)+++-- like <*>+infixl 4 <*>.+, <*>.-, <*>.-$++{- |+> addPair :: (Additive.C a, Additive.C b) => (a,b) -> (a,b) -> (a,b)+> addPair = Elem.run2 $ Elem.with (,) <*>.+  fst <*>.+  snd+-}+(<*>.+) ::+   (C x) =>+   Elem.T (v,v) (x -> a) -> (v -> x) -> Elem.T (v,v) a+(<*>.+) f acc =+   f <*> elementAdd acc++(<*>.-) ::+   (C x) =>+   Elem.T (v,v) (x -> a) -> (v -> x) -> Elem.T (v,v) a+(<*>.-) f acc =+   f <*> elementSub acc++(<*>.-$) ::+   (C x) =>+   Elem.T v (x -> a) -> (v -> x) -> Elem.T v a+(<*>.-$) f acc =+   f <*> elementNeg acc+++-- * Instances for atomic types+ instance C Integer where    {-# INLINE zero #-}    {-# INLINE negate #-}@@ -224,27 +286,27 @@   -{-* Instances for composed types -}+-- * Instances for composed types  instance (C v0, C v1) => C (v0, v1) where    {-# INLINE zero #-}    {-# INLINE negate #-}    {-# INLINE (+) #-}    {-# INLINE (-) #-}-   zero                   = (zero, zero)-   (+)    (x0,x1) (y0,y1) = ((+) x0 y0, (+) x1 y1)-   (-)    (x0,x1) (y0,y1) = ((-) x0 y0, (-) x1 y1)-   negate (x0,x1)         = (negate x0, negate x1)+   zero   = (,) zero zero+   (+)    = Elem.run2 $ pure (,) <*>.+  fst <*>.+  snd+   (-)    = Elem.run2 $ pure (,) <*>.-  fst <*>.-  snd+   negate = Elem.run  $ pure (,) <*>.-$ fst <*>.-$ snd  instance (C v0, C v1, C v2) => C (v0, v1, v2) where    {-# INLINE zero #-}    {-# INLINE negate #-}    {-# INLINE (+) #-}    {-# INLINE (-) #-}-   zero                         = (zero, zero, zero)-   (+)    (x0,x1,x2) (y0,y1,y2) = ((+) x0 y0, (+) x1 y1, (+) x2 y2)-   (-)    (x0,x1,x2) (y0,y1,y2) = ((-) x0 y0, (-) x1 y1, (-) x2 y2)-   negate (x0,x1,x2)            = (negate x0, negate x1, negate x2)+   zero   = (,,) zero zero zero+   (+)    = Elem.run2 $ pure (,,) <*>.+  fst3 <*>.+  snd3 <*>.+  thd3+   (-)    = Elem.run2 $ pure (,,) <*>.-  fst3 <*>.-  snd3 <*>.-  thd3+   negate = Elem.run  $ pure (,,) <*>.-$ fst3 <*>.-$ snd3 <*>.-$ thd3   instance (C v) => C [v] where@@ -268,7 +330,7 @@    (-)    f g x = (-) (f x) (g x)    negate f   x = negate (f x) -{- * Properties -}+-- * Properties  propAssociative :: (Eq a, C a) => a -> a -> a -> Bool propCommutative :: (Eq a, C a) => a -> a -> Bool
+ src/Algebra/AffineSpace.hs view
@@ -0,0 +1,247 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-}+{- |+This module is not yet exported+since its interface is not mature.+There are two approaches for representing affine spaces:++[1] Two sets: A set of points and a set of vectors.+    Examples: Absolute potential and voltage,+    absolute temperature and temperature difference.+    Operations are+      add :: vector -> point -> point+      sub :: point -> point -> vector++[2] One set for points, no vectors.+    Examples: Interpolation+    Operation:+      combine :: [(coefficient, vector)] -> vector+    Where it must be asserted,+    that the coefficients sum up to 1.++The second one is the one we follow here.+It is more similar to Module and VectorSpace.+-}+module Algebra.AffineSpace where++import qualified Algebra.PrincipalIdealDomain as PID+import qualified Algebra.Additive as Additive+import qualified Algebra.Module as Module+import qualified Number.Ratio as Ratio++import qualified Number.Complex as Complex++import Control.Applicative (Applicative(pure, (<*>)), )++import NumericPrelude hiding (zero, )+import PreludeBase+import Prelude ()++{- |+The type class is for representing affine spaces via affine combinations.+However, we didn't find a way to both ensure the property+that the combination coefficients sum up to 1,+and keep it efficient.++We propose this class instead of a combination of Additive and Module+for interpolation for those types,+where scaling and addition alone makes no sense.+Such types are e.g. internal filter parameters in signal processing:+For these types interpolation makes definitely sense,+but addition and scaling make not.++That is, both classes are isomorphic+(you can define one in terms of the other),+but instances of this class are more easily defined,+and using an AffineSpace constraint instead of Module in a type signature+is important for documentation purposes.+AffineSpace should be superclass of Module.+(But then you may ask, why not adding another superclass for Convex spaces.+This class would provide a linear combination operation,+where the coefficients sum up to one+and all of them are non-negative.)++We may add a safety layer that ensures+that the coefficients sum up to 1,+using start points on the simplex+and functions to move on the simplex.+Start points have components that sum up to 1, e.g.++> (1, 0, 0, 0)+> (0, 1, 0, 0)+> (0, 0, 1, 0)+> (0, 0, 0, 1)+> (1/4, 1/4, 1/4, 1/4)++Then you may move along the simplex in the directions++> (1,  -1, 0,  0)+> (0,   1, 0, -1)+> (-1, -1, 3, -1)++which are characterized by components that sum up to 0.++For example linear combination is defined by++> lerp k (a,b) = (1-k)*>a + k*>b++that is the coefficients are (1-k) and k.+The pair (1-k, k) can be constructed+by starting at pair (1,0)+and moving k units in direction (-1,1).++> (1-k, k) = (1,0) + k*(-1,1)++It is however a challenge to manage the coefficient tuples+in a type safe and efficient way.+For small numbers of interpolation nodes+(constant, linear, cubic interpolation)+a type level list would appropriate,+but what to do for large tuples+like for Whittaker interpolation?+++As side note:+In principle it would be sufficient+to provide an affine combination of two points,+since all affine combinations of more points+can be decomposed into such simple combinations.++> lerp a x y = (1-a)*>x + a*>y++E.g. @a*>x + b*>y + c*>z@ with @a+b+c=1@+can be written as @lerp c (lerp (b/(1-c)) x y) z@.+More generally you can use++> lerpnorm a b x y+>    = lerp (b/(a+b)) x y+>    = (a/(a+b))*>x + (b/(a+b))*>y++for writing++> a*>x + b*>y + c*>z ==+>    lerpnorm (a+b) c (lerpnorm a b x y) z++or++> a*>x + b*>y + c*>z + d*>w ==+>    lerpnorm (a+b+c) d (lerpnorm (a+b) c (lerpnorm a b x y) z) w++with @a+b+c+d=1@.++The downside is, that lerpnorm requires division, that is, a field,+whereas the computation of the coefficients+sometimes only requires ring operations.+-}+class Zero v => C a v where+   multiplyAccumulate :: (a,v) -> v -> v++class Zero v where+   zero :: v+++instance Zero Float where+   {-# INLINE zero #-}+   zero = Additive.zero++instance Zero Double where+   {-# INLINE zero #-}+   zero = Additive.zero++instance (Zero a) => Zero (Complex.T a) where+   {-# INLINE zero #-}+   zero = zero Complex.+: zero++instance (PID.C a) => Zero (Ratio.T a) where+   {-# INLINE zero #-}+   zero = Additive.zero+++instance C Float Float where+   {-# INLINE multiplyAccumulate #-}+   multiplyAccumulate (a,x) y = a*x+y++instance C Double Double where+   {-# INLINE multiplyAccumulate #-}+   multiplyAccumulate (a,x) y = a*x+y++instance (C a v) => C a (Complex.T v) where+   {-# INLINE multiplyAccumulate #-}+   multiplyAccumulate =+      makeMac2 (Complex.+:) Complex.real Complex.imag++instance (PID.C a) => C (Ratio.T a) (Ratio.T a) where+   {-# INLINE multiplyAccumulate #-}+   multiplyAccumulate (a,x) y = a*x+y+++infixl 6 *.+++{- |+Infix variant of 'multiplyAccumulate'.+-}+{-# INLINE (*.+) #-}+(*.+) :: C a v => v -> (a,v) -> v+(*.+) = flip multiplyAccumulate+++-- * convenience functions for defining multiplyAccumulate++{-# INLINE multiplyAccumulateModule #-}+multiplyAccumulateModule ::+   Module.C a v =>+   (a,v) -> v -> v+multiplyAccumulateModule (a,x) y =+   a *> x + y+++{- |+A special reader monad.+-}+newtype MAC a v x = MAC {runMac :: (a,v) -> v -> x}++{-# INLINE element #-}+element ::+   (C a x) =>+   (v -> x) -> MAC a v x+element f =+   MAC (\(a,x) y -> multiplyAccumulate (a, f x) (f y))++instance Functor (MAC a v) where+   {-# INLINE fmap #-}+   fmap f (MAC x) =+      MAC $ \av v -> f $ x av v++instance Applicative (MAC a v) where+   {-# INLINE pure #-}+   {-# INLINE (<*>) #-}+   pure x = MAC $ \ _av _v -> x+   MAC f <*> MAC x =+      MAC $ \av v -> f av v $ x av v++{-# INLINE makeMac #-}+makeMac ::+   (C a x) =>+   (x -> v) ->+   (v -> x) ->+   (a,v) -> v -> v+makeMac cons x =+   runMac $ pure cons <*> element x++{-# INLINE makeMac2 #-}+makeMac2 ::+   (C a x, C a y) =>+   (x -> y -> v) ->+   (v -> x) -> (v -> y) ->+   (a,v) -> v -> v+makeMac2 cons x y =+   runMac $ pure cons <*> element x <*> element y++{-# INLINE makeMac3 #-}+makeMac3 ::+   (C a x, C a y, C a z) =>+   (x -> y -> z -> v) ->+   (v -> x) -> (v -> y) -> (v -> z) ->+   (a,v) -> v -> v+makeMac3 cons x y z =+   runMac $ pure cons <*> element x <*> element y <*> element z
src/Algebra/Module.hs view
@@ -23,23 +23,31 @@ import qualified Algebra.Laws as Laws  import Algebra.Ring     ((*), fromInteger, )-import Algebra.Additive ((+), zero, )+import Algebra.Additive ((+), zero, sum, ) +import qualified NumericPrelude.Elementwise as Elem+import Control.Applicative (Applicative(pure, (<*>)), )+ import Data.Function.HT (powerAssociative, )-import Data.List (map, zipWith, foldl, )+import Data.List (map, zipWith, )+import Data.Tuple.HT (fst3, snd3, thd3, )+import Data.Tuple (fst, snd, ) -import Prelude((.), Eq, Bool, Int, Integer, Float, Double)+import Prelude((.), Eq, Bool, Int, Integer, Float, Double, ($), ) -- import qualified Prelude as P   -- Is this right? infixr 7 *> -{- Functional dependency can't be used-   since the instance (Algebra.Module.C a a)-   would conflict with all others.-   class Algebra.Module.C b a | b -> a where -}+{-+Functional dependency can't be used+since @Complex.T a@ is a module+with respect to both @a@ and @Complex.T a@. +class Algebra.Module.C a v | v -> a where+-}+ {-| A Module over a ring satisfies: @@ -47,10 +55,19 @@ >   (a * b) *> c === a *> (b *> c) >   (a + b) *> c === a *> c + b *> c -}-class (Additive.C b, Ring.C a) => C a b where+class (Ring.C a, Additive.C v) => C a v where     -- | scale a vector by a scalar-    (*>) :: a -> b -> b+    (*>) :: a -> v -> v ++(<*>.*>) ::+   (C a x) =>+   Elem.T (a,v) (x -> c) -> (v -> x) -> Elem.T (a,v) c+(<*>.*>) f acc =+   f <*> Elem.element (\(a,v) -> a *> acc v)+++ {-* Instances for atomic types -}  instance C Float Float where@@ -83,17 +100,19 @@  instance (C a b0, C a b1) => C a (b0, b1) where    {-# INLINE (*>) #-}-   s *> (x0,x1)   = (s *> x0, s *> x1)+   (*>) = Elem.run2 $ pure (,) <*>.*> fst <*>.*> snd+   -- s *> (x0,x1)   = (s *> x0, s *> x1)  instance (C a b0, C a b1, C a b2) => C a (b0, b1, b2) where    {-# INLINE (*>) #-}-   s *> (x0,x1,x2) = (s *> x0, s *> x1, s *> x2)+   (*>) = Elem.run2 $ pure (,,) <*>.*> fst3 <*>.*> snd3 <*>.*> thd3+   -- s *> (x0,x1,x2) = (s *> x0, s *> x1, s *> x2) -instance (C a b) => C a [b] where+instance (C a v) => C a [v] where    {-# INLINE (*>) #-}    (*>) = map . (*>) -instance (C a b) => C a (c -> b) where+instance (C a v) => C a (c -> v) where    {-# INLINE (*>) #-}    (*>) s f = (*>) s . f @@ -106,8 +125,8 @@ ToDo: Should it use 'NumericPrelude.List.zipWithMatch' ? -}-linearComb :: C a b => [a] -> [b] -> b-linearComb c = foldl (+) zero . zipWith (*>) c+linearComb :: C a v => [a] -> [v] -> v+linearComb c = sum . zipWith (*>) c  {-| This function can be used to define any@@ -116,18 +135,18 @@ Better move to "Algebra.Additive"? -} {-# INLINE integerMultiply #-}-integerMultiply :: (ToInteger.C a, Additive.C b) => a -> b -> b-integerMultiply a b =-   powerAssociative (+) zero b (ToInteger.toInteger a)+integerMultiply :: (ToInteger.C a, Additive.C v) => a -> v -> v+integerMultiply a v =+   powerAssociative (+) zero v (ToInteger.toInteger a)   {- * Properties -} -propCascade :: (Eq b, C a b) => b -> a -> a -> Bool+propCascade :: (Eq v, C a v) => v -> a -> a -> Bool propCascade  =  Laws.leftCascade (*) (*>) -propRightDistributive :: (Eq b, C a b) => a -> b -> b -> Bool+propRightDistributive :: (Eq v, C a v) => a -> v -> v -> Bool propRightDistributive  =  Laws.rightDistributive (*>) (+) -propLeftDistributive :: (Eq b, C a b) => b -> a -> a -> Bool+propLeftDistributive :: (Eq v, C a v) => v -> a -> a -> Bool propLeftDistributive x  =  Laws.homomorphism (*>x) (+) (+)
src/MathObj/Gaussian/Bell.hs view
@@ -20,40 +20,50 @@ import Algebra.Additive ((+), )  import Test.QuickCheck (Arbitrary, arbitrary, coarbitrary, )-import Control.Monad (liftM3, )+import Control.Monad (liftM4, )  -- import Prelude (($)) import NumericPrelude import PreludeBase hiding (reverse, )  -data T a = Cons {c0, c1 :: Complex.T a, c2 :: a}+data T a = Cons {amp :: a, c0, c1 :: Complex.T a, c2 :: a}    deriving (Eq, Show)  instance (Real.C a, Arbitrary a) => Arbitrary (T a) where    arbitrary =-      liftM3-         (\a b c -> Cons a b (1 + abs c))-         arbitrary arbitrary arbitrary+      liftM4+         (\k a b c -> Cons k a b (1 + abs c))+         arbitrary arbitrary arbitrary arbitrary    coarbitrary = undefined  -constant :: Additive.C a => T a-constant = Cons zero zero zero+constant :: Ring.C a => T a+constant = Cons one zero zero zero +{- |+eigenfunction of 'fourier'+-}+unit :: Ring.C a => T a+unit = Cons one zero zero one+ {-# INLINE evaluate #-} evaluate :: (Trans.C a) =>    T a -> a -> Complex.T a evaluate f x =-   Complex.exp $ Complex.scale (-pi) $-   c0 f + Complex.scale x (c1 f) + Complex.fromReal (c2 f * x^2)+   Complex.scale+     (sqrt (amp f))+     (Complex.exp $ Complex.scale (-pi) $+      c0 f + Complex.scale x (c1 f) + Complex.fromReal (c2 f * x^2))  evaluateSqRt :: (Trans.C a) =>    T a -> a -> Complex.T a evaluateSqRt f x0 =-   let x = sqrt pi * x0-   in  Complex.exp $ negate $-       c0 f + Complex.scale x (c1 f) + Complex.fromReal (c2 f * x^2)+   Complex.scale+     (sqrt (amp f))+     (let x = sqrt pi * x0+      in  Complex.exp $ negate $+          c0 f + Complex.scale x (c1 f) + Complex.fromReal (c2 f * x^2))  exponentPolynomial :: (Additive.C a) =>    T a -> Poly.T (Complex.T a)@@ -61,10 +71,17 @@    Poly.fromCoeffs [c0 f, c1 f, Complex.fromReal (c2 f)]  -multiply :: (Additive.C a) =>+variance :: (Trans.C a) =>+   T a -> a+variance f =+   recip $ c2 f * 2*pi++multiply :: (Ring.C a) =>    T a -> T a -> T a multiply f g =-   Cons (c0 f + c0 g) (c1 f + c1 g) (c2 f + c2 g)+   Cons+      (amp f * amp g)+      (c0 f + c0 g) (c1 f + c1 g) (c2 f + c2 g)   {-@@ -89,6 +106,7 @@          = -(f1 - g1)^2/(4*(f2 + g2))        -}    in  Cons+          ((amp f * amp g) / (c2 f + c2 g))           (c0 f + c0 g               - Complex.scale (recip (4*s)) ((c1 f - c1 g)^2))           (Complex.scale (c2 g / s) (c1 f) +@@ -110,6 +128,7 @@        g1 = translateComplex gd g0    in  translateComplex (negate $ fd + gd) $        Cons+          ((amp f0 * amp g0) / (c2 f0 + c2 g0))           (c0 f1 + c0 g1) zero           (recip $ recip (c2 f1) + recip (c2 g1)) @@ -126,6 +145,7 @@        c = c2 f        rc = recip c    in  Cons+          (amp f / c2 f)           (Complex.scale (rc/4) (-b^2) + a)           (Complex.scale rc $ Complex.quarterRight b)           rc@@ -134,7 +154,7 @@    T a -> T a fourierByTranslation f =    translateComplex (Complex.scale (1/2) $ Complex.quarterLeft $ c1 f) $-   Cons (c0 f) zero (recip $ c2 f)+   Cons (amp f / c2 f) (c0 f) zero (recip $ c2 f)  {- a + b*x + c*x^2@@ -201,6 +221,7 @@        b = c1 f        c = c2 f    in  Cons+          (amp f)           (Complex.fromReal (c*d^2) - Complex.scale d b + a)           (Complex.fromReal (-2*c*d) + b)           c@@ -211,6 +232,7 @@        b = c1 f        c = c2 f    in  Cons+          (amp f)           (Complex.scale c (d^2) - b*d + a)           (Complex.scale (-2*c) d + b)           c@@ -218,6 +240,7 @@ modulate :: Ring.C a => a -> T a -> T a modulate d f =    Cons+      (amp f)       (c0 f)       (c1 f + (zero +: 2*d))       (c2 f)@@ -225,6 +248,7 @@ turn :: Ring.C a => a -> T a -> T a turn d f =    Cons+      (amp f)       (c0 f + (zero +: 2*d))       (c1 f)       (c2 f)@@ -237,6 +261,7 @@ dilate :: Field.C a => a -> T a -> T a dilate k f =    Cons+      (amp f)       (c0 f)       (Complex.scale (recip k) $ c1 f)       (c2 f / k^2)@@ -244,9 +269,18 @@ shrink :: Ring.C a => a -> T a -> T a shrink k f =    Cons+      (amp f)       (c0 f)       (Complex.scale k $ c1 f)       (k^2 * c2 f)++amplify :: (Ring.C a) => a -> T a -> T a+amplify k f =+   Cons+      (k^2 * amp f)+      (c0 f)+      (c1 f)+      (c2 f)   {- laws
+ src/MathObj/Gaussian/Example.hs view
@@ -0,0 +1,227 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-+Reciprocal of variance of a Gaussian bell curve.+We describe the curve only in terms of its variance+thus we represent a bell curve at the coordinate origin+neglecting its amplitude.++We could also define the amplitude as @root 4 c@,+thus preserving L2 norm being one,+but then @dilate@ and @shrink@ also include an amplification.++We could do some projective geometry in the exponent+in order to also have zero variance,+which corresponds to the dirac impulse.+-}+module MathObj.Gaussian.Example where++import qualified MathObj.Gaussian.Polynomial as PolyBell+import qualified MathObj.Gaussian.Bell as Bell+import qualified MathObj.Gaussian.Variance as Var++import qualified MathObj.Polynomial as Poly++import qualified Algebra.Transcendental as Trans+import qualified Algebra.Algebraic      as Algebraic+import qualified Algebra.Field          as Field+-- import qualified Algebra.Real           as Real+import qualified Algebra.Ring           as Ring+-- import qualified Algebra.Additive       as Additive++import qualified Number.Complex as Complex++import Algebra.Transcendental (pi, )+import Algebra.Algebraic (root, )+import Algebra.Ring ((*), (^), )++import Number.Complex ((+:), )++import qualified Numerics.Function as Func+import qualified Numerics.Fourier as Fourier+import qualified Numerics.Integration as Integ+import qualified Numerics.Differentiation as Diff++import qualified Graphics.Gnuplot.Simple as GP++import Control.Applicative (liftA2, )++-- import System.Exit (ExitCode, )++-- import Prelude (($))+import NumericPrelude+import PreludeBase+import qualified Prelude as P+++curve0 :: Var.T Double+curve0 = curve0a++curve0a :: Var.T Double+curve0a = Var.Cons 1.4 3.3++curve0b :: Var.T Double+curve0b = Var.Cons 2.2 1.7++variance0 :: (Double, Double)+variance0 =+   (Var.variance curve0,+    (Integ.rectangular 1000 (-2,2) $ liftA2 (*) (^2) (Var.evaluate curve0)) /+    (Integ.rectangular 1000 (-2,2) $ Var.evaluate curve0))++norm10 :: (Double, Double)+norm10 =+   (Integ.rectangular 1000 (-2,2) $ Var.evaluate curve0,+    Var.norm1 curve0)++norm20 :: (Double, Double)+norm20 =+   (sqrt $ Integ.rectangular 1000 (-2,2) $ (^2) . Var.evaluate curve0,+    Var.norm2 curve0)++norm30 :: (Double, Double)+norm30 =+   (root 3 $ Integ.rectangular 1000 (-2,2) $ (^3) . Var.evaluate curve0,+    Var.normP 3 curve0)++fourier0 :: IO ()+fourier0 =+   GP.plotFuncs []+      (GP.linearScale 100 (-2,2))+      [Var.evaluate $ Var.fourier curve0,+       Fourier.analysisTransformOneReal 100 (-2,2) $ Var.evaluate curve0]++multiply0 :: IO ()+multiply0 =+   GP.plotFuncs []+      (GP.linearScale 100 (-1,1))+      [Var.evaluate $ Var.multiply curve0a curve0b,+       liftA2 (*) (Var.evaluate curve0a) (Var.evaluate curve0b)]++convolve0 :: IO ()+convolve0 =+   GP.plotFuncs []+      (GP.linearScale 100 (-2,2))+      [Var.evaluate $ Var.convolve curve0a curve0b,+       Integ.convolve 1000 (-3,3) (Var.evaluate curve0a) (Var.evaluate curve0b)]+++curve1 :: Bell.T Double+curve1 = curve1a++curve1a :: Bell.T Double+curve1a = Bell.Cons 1.4 (0.1+:0.3) ((-0.2)+:1.4) 2.3++curve1b :: Bell.T Double+curve1b = Bell.Cons 2.2 ((-0.3)+:2.1) (0.2+:(-0.4)) 1.7++variance1 :: (Double, Double)+variance1 =+   (Bell.variance curve1,+    (Integ.rectangular 1000 (-2,2) $+        liftA2 (*) (^2)+           (Complex.magnitudeSqr .+            Func.translateRight+               (Complex.real (Bell.c1 curve1) / (2 * Bell.c2 curve1))+               (Bell.evaluate curve1))) /+    (Integ.rectangular 1000 (-2,2) $ Complex.magnitude . Bell.evaluate curve1))++{- the norm depends on too much things+norm0vs1 :: (Double, Double)+norm0vs1 =+   ((Integ.rectangular 1000 (-5,5) $ Var.evaluate curve0)+         * exp (- Complex.real (Bell.c0 curve1)),+    Integ.rectangular 1000 (-5,5) $ Complex.magnitude . Bell.evaluate curve1)+-}++fourier1 :: IO ()+fourier1 =+   GP.plotFuncs []+      (GP.linearScale 100 (-5,5))+      [Complex.real . (Bell.evaluate $ Bell.fourier curve1),+       fourierAnalysisReal 100 (-2,2) $ Bell.evaluate curve1]+++curve2 :: PolyBell.T Double+curve2 =+   PolyBell.Cons+--      Bell.unit+--      (Bell.Cons 1.4 (0.1+:0.3) 0 1.2)+--      (Bell.Cons 1.4 (0.1+:0.3) ((-0.2)+:1.4) 1)+      curve1+--      (Poly.fromCoeffs [one])+--      (Poly.fromCoeffs [zero,one])+--      (Poly.fromCoeffs [zero,zero,one])+--      (Poly.fromCoeffs [0,Complex.imaginaryUnit])+      (Poly.fromCoeffs [1.4+:(-0.1),0.8+:(0.1),(-1.1)+:0.3])++differentiate2 :: IO ()+differentiate2 =+   GP.plotFuncs []+      (GP.linearScale 100 (-2,2))+      [Complex.real . (PolyBell.evaluateSqRt $ PolyBell.differentiate curve2),+       ((/ sqrt pi) . ) $ Diff.diff (1e-5) $ Complex.real . PolyBell.evaluateSqRt curve2]++fourier2 :: IO ()+fourier2 =+   GP.plotFuncs []+      (GP.linearScale 100 (-5,5))+      [Complex.real . (PolyBell.evaluateSqRt $ PolyBell.fourier curve2),+       fourierAnalysisReal 100 (-2,2) $ PolyBell.evaluateSqRt curve2]++++fourierAnalysisReal ::+   (P.Floating a) =>+   Integer -> (a, a) -> (a -> Complex.T a) -> a -> a+fourierAnalysisReal n rng f =+   liftA2 (P.-)+      (Fourier.analysisTransformOneReal n rng (Complex.real . f))+      (Fourier.analysisTransformOneImag n rng (Complex.imag . f))+++{- |+Try to approximate @\x -> exp (-x^2) * x@+by a difference of translated Gaussian bells.++exp(-x^2) * x+  ==  exp(-(a+b*x+c*x^2)) - exp(-(a-b*x+c*x^2))+  ==  exp(-(a+c*x^2)) * (exp(-b*x) - exp(b*x))+  ==  exp(-(a+c*x^2)) * 2*sinh (b*x)++It holds+  lim (\b x -> sinh (b*x) / b)  =  id+-}+diffApprox :: IO ()+diffApprox =+   let amp = (2*b)^- (-2)+       a = 0+       {-+       amp = 1+       a = log (2 * abs b)+       -}+       b = -0.1+       c = 1+       ac = Complex.fromReal a+       bc = Complex.fromReal b+   in  GP.plotFuncs []+          (GP.linearScale 100 (-2,2::Double))+          [Complex.real .+           (PolyBell.evaluateSqRt $+              PolyBell.Cons Bell.unit (Poly.fromCoeffs [zero,one])),+           Complex.real .+           liftA2 (-)+             (PolyBell.evaluateSqRt $+                PolyBell.Cons (Bell.Cons amp ac bc c) (Poly.fromCoeffs [one]))+             (PolyBell.evaluateSqRt $+                PolyBell.Cons (Bell.Cons amp ac (-bc) c) (Poly.fromCoeffs [one]))]+++polyApprox :: IO ()+polyApprox =+   GP.plotFuncs []+      (GP.linearScale 100 (-2,2::Double))+      [Complex.real .+         PolyBell.evaluateSqRt curve2,+       Complex.real . sum .+         mapM (\(amp,b) -> \x -> amp * Bell.evaluateSqRt b x)+         (PolyBell.approximateByBells 0.1 curve2)]
src/MathObj/Gaussian/Polynomial.hs view
@@ -7,17 +7,30 @@ as unit for translations and modulations, for linear factors and in the differentiation. -}+{-+ToDo:++* In order to avoid the weird @sqrt pi@ factor,+  use a polynomial expression in @pi@.++* sum of multiple bells using Data.Map from exponent polynomial to coefficient polynomial+  use of Algebra object.++* Projective geometry in order to support Dirac impulse.+-} module MathObj.Gaussian.Polynomial where  import qualified MathObj.Gaussian.Bell as Bell +import qualified MathObj.LaurentPolynomial as LPoly import qualified MathObj.Polynomial as Poly import qualified Number.Complex     as Complex +import qualified Algebra.ZeroTestable   as ZeroTestable import qualified Algebra.Differential   as Differential import qualified Algebra.Transcendental as Trans import qualified Algebra.Field          as Field--- import qualified Algebra.Real           as Real+import qualified Algebra.Real           as Real import qualified Algebra.Ring           as Ring import qualified Algebra.Additive       as Additive @@ -25,65 +38,165 @@ import Algebra.Ring ((*), ) -- import Algebra.Additive ((+)) +import qualified Data.Record.HT as Rec+import qualified Data.List as List+import Data.Function.HT (nest, )+import Data.Eq.HT (equating, )+import Data.List.HT (mapAdjacent, )+import Data.Tuple.HT (forcePair, )++import Test.QuickCheck (Arbitrary, arbitrary, coarbitrary, )+import Control.Monad (liftM2, )+ import NumericPrelude import PreludeBase hiding (reverse, ) -- import Prelude ()   data T a = Cons {bell :: Bell.T a, polynomial :: Poly.T (Complex.T a)}-   deriving (Eq, Show)+   deriving (Show) +instance Real.C a => Eq (T a) where+   (==) = equal -{-# INLINE evaluate #-}-evaluate :: (Trans.C a) =>+{-+The derived Eq is not correct.+We have to combine the amplitude of the bell with the polynomial,+respecting signs and the square root of the bell amplitude.+-}+equal :: Real.C a => T a -> T a -> Bool+equal x y =+   let bx = bell x+       by = bell y+       csign c =+          Complex.real c > 0 ||+          (Complex.real c == 0 && Complex.imag c > 0)+       scaleSqr b =+          map (\c -> (Complex.scale (Bell.amp b) (c^2), csign c)) .+          Poly.coeffs . polynomial+   in  Rec.equal+          (equating Bell.c0 :+           equating Bell.c1 :+           equating Bell.c2 :+           [])+          bx by+       &&+       scaleSqr by x == scaleSqr bx y+++instance (Real.C a, Arbitrary a) => Arbitrary (T a) where+   arbitrary =+--      liftM2 Cons arbitrary arbitrary+      liftM2 Cons+         arbitrary+         -- we have to restrict the number of polynomial coefficients,+         -- since with the quadratic time algorithms like fourier and convolve,+         -- in connection with Rational slow down tests too much.+         (fmap (Poly.fromCoeffs . take 5 . Poly.coeffs) arbitrary)+   coarbitrary = undefined++++{-# INLINE evaluateSqRt #-}+evaluateSqRt :: (Trans.C a) =>    T a -> a -> Complex.T a-evaluate f x =+evaluateSqRt f x =    Bell.evaluateSqRt (bell f) x *    Poly.evaluate (polynomial f) (Complex.fromReal $ sqrt pi * x) {- ToDo: evaluating a complex polynomial for a real argument can be optimized -}  +constant :: (Ring.C a) => T a+constant =+   Cons Bell.constant (Poly.const one)++scale :: (Ring.C a) => a -> T a -> T a+scale x f =+   f{polynomial = fmap (Complex.scale x) $ polynomial f}++scaleComplex :: (Ring.C a) => Complex.T a -> T a -> T a+scaleComplex x f =+   f{polynomial = fmap (x*) $ polynomial f}+++eigenfunction :: (Field.C a) => Int -> T a+eigenfunction =+   eigenfunctionDifferential++eigenfunction0 :: (Ring.C a) => T a+eigenfunction0 =+   Cons Bell.unit (Poly.fromCoeffs [one])++eigenfunction1 :: (Ring.C a) => T a+eigenfunction1 =+   Cons Bell.unit (Poly.fromCoeffs [zero, one])++eigenfunction2 :: (Field.C a) => T a+eigenfunction2 =+   Cons Bell.unit (Poly.fromCoeffs [-(1/4), zero, one])++eigenfunction3 :: (Field.C a) => T a+eigenfunction3 =+   Cons Bell.unit (Poly.fromCoeffs [zero, -(3/4), zero, one])+++eigenfunctionDifferential :: (Field.C a) => Int -> T a+eigenfunctionDifferential n =+   (\f -> f{bell = Bell.unit}) $+   nest n (scale (-1/4) . differentiate) $+   Cons (Bell.Cons one zero zero 2) one++eigenfunctionIterative :: (Field.C a, Real.C a) => Int -> T a+eigenfunctionIterative n =+   fst . head . dropWhile (uncurry (/=)) . mapAdjacent (,) $+   eigenfunctionIteration $+   Cons+      Bell.unit+      (Poly.fromCoeffs $ replicate n zero ++ [one])++eigenfunctionIteration :: (Field.C a) => T a -> [T a]+eigenfunctionIteration =+   iterate (\x ->+      let y = fourier x+          px = polynomial x+          py = polynomial y+          c = last (Poly.coeffs px) / last (Poly.coeffs py)+      in  y{polynomial = fmap (0.5*) (px + fmap (c*) py)})++ multiply :: (Ring.C a) =>    T a -> T a -> T a-multiply x y =+multiply f g =    Cons-      (Bell.multiply (bell x) (bell y))-      (polynomial x * polynomial y)+      (Bell.multiply (bell f) (bell g))+      (polynomial f * polynomial g)  convolve :: (Field.C a) =>    T a -> T a -> T a convolve f g =    reverse $ fourier $ multiply (fourier f) (fourier g) -reverse :: Additive.C a => T a -> T a-reverse x =-   Cons-      (Bell.reverse $ bell x)-      (Poly.reverse $ polynomial x)- {- We use a Horner like scheme in order to translate multiplications with @id@ to differentations on the Fourier side. Quadratic runtime. -fourier (Cons bell (Poly.const a + Poly.shift x))-  = fourier (Cons bell (Poly.const a)) + fourier (Cons bell (Poly.shift x))-  = fourier (Cons bell (Poly.const a)) + differentiate (fourier (Cons bell x))--untested+fourier (Cons bell (Poly.const a + Poly.shift f))+  = fourier (Cons bell (Poly.const a)) + fourier (Cons bell (Poly.shift f))+  = fourier (Cons bell (Poly.const a)) + C * differentiate (fourier (Cons bell f)) -} fourier :: (Field.C a) =>    T a -> T a-fourier x =+fourier f =    foldr       (\c p ->           let q = differentiate p           in  q{polynomial =                    Poly.const c +-                   fmap Complex.quarterLeft (polynomial q)})-      (Cons (Bell.fourier $ bell x) zero) $-   Poly.coeffs $ polynomial x+                   fmap (Complex.scale (1/2) . Complex.quarterLeft) (polynomial q)})+      (Cons (Bell.fourier $ bell f) zero) $+   Poly.coeffs $ polynomial f  {- Differentiate and divide by @sqrt pi@ in order to stay in a ring.@@ -92,8 +205,174 @@ differentiate :: (Ring.C a) => T a -> T a differentiate f =    f{polynomial =-        Bell.exponentPolynomial (bell f) * polynomial f +-        Differential.differentiate (polynomial f)}+        Differential.differentiate (polynomial f)+        - Differential.differentiate (Bell.exponentPolynomial (bell f))+           * polynomial f}++{-+snd $ integrate $ differentiate (Cons Bell.unit (Poly.fromCoeffs [7,7,7,7]) :: T Double)+-}+integrate ::+   (Field.C a, ZeroTestable.C a) =>+   T a -> (Complex.T a, T a)+integrate f =+   let fs = Poly.coeffs $ polynomial f+       (ys,~[r]) =+          Poly.divModRev+             {-+             We need the shortening convention of 'zipWith'+             in order to limit the result list,+             we cannot use list instance for (-).+             -}+             (zipWith (-)+                (0 : 0 : diffRev ys)+                (List.reverse fs))+             (List.reverse $ Poly.coeffs $+              Differential.differentiate $+              Bell.exponentPolynomial $ bell f)+   in  forcePair $+       if null fs+         then (zero, f)+         else (r, f{polynomial = Poly.fromCoeffs $ List.reverse ys})++diffRev :: Ring.C a => [a] -> [a]+diffRev xs =+   zipWith (*) xs+      (drop 1 (iterate (subtract 1) (fromIntegral $ length xs)))++translate :: Ring.C a => a -> T a -> T a+translate d =+   translateComplex (Complex.fromReal d)++translateComplex :: Ring.C a => Complex.T a -> T a -> T a+translateComplex d f =+   Cons+      (Bell.translateComplex d $ bell f)+      (Poly.translate d $ polynomial f)++modulate :: Ring.C a => a -> T a -> T a+modulate d f =+   Cons+      (Bell.modulate d $ bell f)+      (polynomial f)++turn :: Ring.C a => a -> T a -> T a+turn d f =+   Cons+      (Bell.turn d $ bell f)+      (polynomial f)++reverse :: Additive.C a => T a -> T a+reverse f =+   Cons+      (Bell.reverse $ bell f)+      (Poly.reverse $ polynomial f)++dilate :: Field.C a => a -> T a -> T a+dilate k f =+   Cons+      (Bell.dilate k $ bell f)+      (Poly.dilate (Complex.fromReal k) $ polynomial f)++shrink :: Ring.C a => a -> T a -> T a+shrink k f =+   Cons+      (Bell.shrink k $ bell f)+      (Poly.shrink (Complex.fromReal k) $ polynomial f)++{-+We could also amplify the polynomial coefficients.+-}+amplify :: Ring.C a => a -> T a -> T a+amplify k f =+   Cons+      (Bell.amplify k $ bell f)+      (polynomial f)+++{- |+Approximate a @T a@ using a linear combination of translated @Bell.T a@.+The smaller the unit (e.g. 0.1, 0.01, 0.001)+the better the approximation but the worse the numeric properties.++We cannot put all information into @amp@ of @Bell@,+since @amp@ must be real, but is complex here by construction.+We really need at least signed amplitudes at this place,+since we want to represent differences of Gaussians.+-}+approximateByBells ::+   Field.C a =>+   a -> T a -> [(Complex.T a, Bell.T a)]+approximateByBells unit f =+   let b = bell f+       amps =+          -- approximateByBellsByTranslation+          approximateByBellsAtOnce+             unit+             (Complex.scale (recip (2 * Bell.c2 b)) (Bell.c1 b))+             (recip (2*unit*Bell.c2 b))+             (polynomial f)+   in  zip (LPoly.coeffs amps) $+       map+          (\d -> Bell.translate d b)+          (laurentAbscissas (unit/2) amps)++approximateByBellsAtOnce ::+   Field.C a =>+   a -> Complex.T a -> a -> Poly.T (Complex.T a) -> LPoly.T (Complex.T a)+approximateByBellsAtOnce unit d s p =+   foldr+      (\x amps0 ->+         {-+         Decompose (bell t * (t-d)) = bell t * t - bell t * d+         -}+         let y = fmap (Complex.scale s) amps0+         in  -- \t -> bell t * t+             --    ~   (translate unit bell - translate (-unit) bell) / unit+             LPoly.shift 1 y -+             LPoly.shift (-1) y ++             -- bell t * d+             zipWithAbscissas+                (\t z -> (Complex.fromReal t - d) * z)+                (unit/2) amps0 ++             LPoly.const x)+      (LPoly.fromCoeffs [])+      (Poly.coeffs p)++approximateByBellsByTranslation ::+   Field.C a =>+   a -> Complex.T a -> a -> Poly.T (Complex.T a) -> LPoly.T (Complex.T a)+approximateByBellsByTranslation unit d s p =+   foldr+      (\x amps0 ->+         {-+         Decompose (bell t * (t-d)) = bell t * t - bell t * d+         -}+         let y = fmap (Complex.scale s) amps0+         in  -- \t -> bell t * t+             --    ~   (translate unit bell - translate (-unit) bell) / unit+             LPoly.shift 1 y -+             LPoly.shift (-1) y ++             -- bell t * d+             zipWithAbscissas Complex.scale (unit/2) amps0 ++             LPoly.const x)+      (LPoly.fromCoeffs [])+      (Poly.coeffs $ Poly.translate d p)++zipWithAbscissas ::+   (Ring.C a) =>+   (a -> b -> c) -> a -> LPoly.T b -> LPoly.T c+zipWithAbscissas h unit y =+   LPoly.fromShiftCoeffs (LPoly.expon y) $+   zipWith h+      (laurentAbscissas unit y)+      (LPoly.coeffs y)++laurentAbscissas :: Ring.C a => a -> LPoly.T c -> [a]+laurentAbscissas unit =+   map (\d -> fromIntegral d * unit) .+   iterate (1+) . LPoly.expon+  {- No Ring instance for Gaussians instance (Ring.C a) => Differential.C (T a) where
src/MathObj/Gaussian/Variance.hs view
@@ -1,12 +1,7 @@ {-# LANGUAGE NoImplicitPrelude #-} {--Reciprocal of variance of a Gaussian bell curve.-We describe the curve only in terms of its variance-thus we represent a bell curve at the coordinate origin-neglecting its amplitude.--We could also define the amplitude as @root 4 c@,-but then @dilate@ and @shrink@ also include an amplification.+We represent a Gaussian bell curve in terms of the reciprocal of its variance+and its value at the origin.  We could do some projective geometry in the exponent in order to also have zero variance,@@ -23,34 +18,38 @@ import qualified Algebra.Ring           as Ring import qualified Algebra.Additive       as Additive +{- import Algebra.Transcendental (pi, ) import Algebra.Ring ((*), (^), ) import Algebra.Additive ((+))-+-} import Test.QuickCheck (Arbitrary, arbitrary, coarbitrary, )-+import Control.Monad (liftM2, )  -- import Prelude (($)) import NumericPrelude import PreludeBase  -data T a = Cons {c :: a}+data T a = Cons {amp, c :: a}    deriving (Eq, Show)  instance (Real.C a, Arbitrary a) => Arbitrary (T a) where-   arbitrary = fmap (Cons . (1+) . abs) arbitrary+   arbitrary =+      liftM2 Cons+         arbitrary+         (fmap ((1+) . abs) arbitrary)    coarbitrary = undefined  -constant :: Additive.C a => T a-constant = Cons zero+constant :: Ring.C a => T a+constant = Cons one zero  {-# INLINE evaluate #-} evaluate :: (Trans.C a) =>    T a -> a -> a evaluate f x =-   exp $ (-pi * c f * x^2)+   sqrt (amp f) * exp (-pi * c f * x^2)  exponentPolynomial :: (Additive.C a) =>    T a -> Poly.T a@@ -58,17 +57,17 @@    Poly.fromCoeffs [zero, zero, c f]  -norm1 :: (Algebraic.C a) => T a -> a+norm1 :: (Algebraic.C a, Real.C a) => T a -> a norm1 f =-   recip $ sqrt $ c f+   sqrt $ abs (amp f) / c f -norm2 :: (Algebraic.C a) => T a -> a+norm2 :: (Algebraic.C a, Real.C a) => T a -> a norm2 f =-   recip $ sqrt $ sqrt $ 2 * c f+   sqrt $ abs (amp f) / (sqrt $ 2 * c f) -normP :: (Trans.C a) => a -> T a -> a+normP :: (Trans.C a, Real.C a) => a -> T a -> a normP p f =-   (p * c f) ^? (- recip (2*p))+   sqrt (abs (amp f)) * (p * c f) ^? (- recip (2*p))   variance :: (Trans.C a) =>@@ -76,10 +75,10 @@ variance f =    recip $ c f * 2*pi -multiply :: (Additive.C a) =>+multiply :: (Ring.C a) =>    T a -> T a -> T a multiply f g =-   Cons $ c f + c g+   Cons (amp f * amp g) (c f + c g)  {- | > convolve x y t =@@ -88,7 +87,9 @@ convolve :: (Field.C a) =>    T a -> T a -> T a convolve f g =-   Cons $ recip $ recip (c f) + recip (c g)+   Cons+      (amp f * amp g / (c f + c g))+      (recip $ recip (c f) + recip (c g))  {- | > fourier x f =@@ -97,18 +98,23 @@ fourier :: (Field.C a) =>    T a -> T a fourier f =-   Cons $ recip $ c f+   Cons (amp f / c f) (recip $ c f) {- fourier (t -> exp(-(a*t)^2)) -}  dilate :: (Field.C a) => a -> T a -> T a dilate k f =-   Cons $ c f / k^2+   Cons (amp f) $ c f / k^2  shrink :: (Ring.C a) => a -> T a -> T a shrink k f =-   Cons $ c f * k^2+   Cons (amp f) $ c f * k^2++amplify :: (Ring.C a) => a -> T a -> T a+amplify k f =+   Cons (k^2 * amp f) $ c f+  {- laws fourier (convolve f g) = multiply (fourier f) (fourier g)
src/MathObj/LaurentPolynomial.hs view
@@ -69,8 +69,12 @@ bounds :: T a -> (Int, Int) bounds (Cons xt x) = (xt, xt + length x - 1) +shift :: Int -> T a -> T a+shift t (Cons xt x) = Cons (xt+t) x++{-# DEPRECATED translate "In order to avoid confusion with Polynomial.translate, use 'shift' instead" #-} translate :: Int -> T a -> T a-translate t (Cons xt x) = Cons (xt+t) x+translate = shift   instance Functor T where
src/MathObj/Polynomial.hs view
@@ -50,11 +50,12 @@     compose, equal, add, sub, negate,     horner, hornerCoeffVector, hornerArgVector,     shift, unShift,-    mul, scale, divMod,+    mul, scale, divMod, divModRev,     tensorProduct, tensorProductAlt,     mulShear, mulShearTranspose,     progression, differentiate, integrate, integrateInt,-    fromRoots, alternate, reverse, )+    fromRoots, alternate, reverse,+    translate, dilate, shrink, ) where  import qualified Algebra.Differential         as Differential@@ -76,8 +77,9 @@ import Control.Monad (liftM, ) import qualified Data.List as List import NumericPrelude.List (zipWithOverlap, )+import Data.Tuple.HT (mapPair, mapFst, forcePair, ) import Data.List.HT-          (dropWhileRev, shear, shearTranspose, outerProduct, )+          (dropWhileRev, switchL, shear, shearTranspose, outerProduct, )  import Test.QuickCheck (Arbitrary(arbitrary,coarbitrary)) @@ -307,19 +309,27 @@  divMod :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a]) divMod x y =-    let (y0:ys) = dropWhile isZero (List.reverse y)-        aux l xs' =-          if l < 0-            then ([], xs')-            else+   mapPair (List.reverse, List.reverse) $+   divModRev (List.reverse x) (List.reverse y)++{-+snd $ Poly.divMod (repeat (1::Double)) [1,1]+-}+divModRev :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])+divModRev x y =+   let (y0:ys) = dropWhile isZero y+       -- the second parameter represents lazily (length x - length y)+       aux xs' =+         forcePair .+         switchL+           ([], xs')+           (P.const $               let (x0:xs) = xs'                   q0      = x0/y0-                  (d',m') = aux (l-1) (sub xs (scale q0 ys))-              in  (q0:d',m')-        (d, m) = aux (length x - length y) (List.reverse x)-    in  if isZero y-          then error "MathObj.Polynomial: division by zero"-          else (List.reverse d, List.reverse m)+              in  mapFst (q0:) . aux (sub xs (scale q0 ys)))+   in  if isZero y+         then error "MathObj.Polynomial: division by zero"+         else aux x (drop (length y - 1) x)  instance (ZeroTestable.C a, Field.C a) => Integral.C (T a) where   divMod (Cons x) (Cons y) =@@ -351,7 +361,7 @@  {-# INLINE differentiate #-} differentiate :: (Ring.C a) => [a] -> [a]-differentiate = zipWith (*) progression . tail+differentiate = zipWith (*) progression . drop 1  {-# INLINE integrate #-} integrate :: (Field.C a) => a -> [a] -> [a]@@ -388,6 +398,17 @@ {-# INLINE reverse #-} reverse :: Additive.C a => T a -> T a reverse = lift1 alternate++translate :: Ring.C a => a -> T a -> T a+translate d =+   lift1 $ foldr (\c p -> [c] + mulLinearFactor d p) []++shrink :: Ring.C a => a -> T a -> T a+shrink k =+   lift1 $ zipWith (*) (iterate (k*) one)++dilate :: Field.C a => a -> T a -> T a+dilate = shrink . Field.recip  {- see htam: Wavelet/DyadicResultant
src/Number/Complex.hs view
@@ -34,6 +34,7 @@         signum,         toPolar,         magnitude,+        magnitudeSqr,         phase,         -- * Conjugate         conjugate,@@ -69,9 +70,16 @@ import qualified Algebra.Indexable          as Indexable  import Algebra.ZeroTestable(isZero)-import Algebra.Module((*>))-import Algebra.Algebraic((^/))+import Algebra.Module((*>), (<*>.*>), )+import Algebra.Algebraic((^/), ) +import qualified NumericPrelude.Elementwise as Elem+import Algebra.Additive ((<*>.+), (<*>.-), (<*>.-$), )++import Foreign.Storable (Storable (..), )+import qualified Foreign.Storable.Record as Store+import Control.Applicative (liftA2, )+ import Test.QuickCheck (Arbitrary, arbitrary, coarbitrary) import Control.Monad (liftM2) @@ -120,7 +128,23 @@    {-# INLINE coarbitrary #-}    coarbitrary = undefined +instance (Storable a) => Storable (T a) where+   sizeOf    = Store.sizeOf store+   alignment = Store.alignment store+   peek      = Store.peek store+   poke      = Store.poke store +store ::+   (Storable a) =>+   Store.Dictionary (T a)+store =+   Store.run $+   liftA2 (+:)+      (Store.element real)+      (Store.element imag)+++ {- * Functions -}  -- | Construct a complex number from real and imaginary part.@@ -265,13 +289,13 @@     {-# SPECULATE instance Additive.C (T Float) #-}     {-# SPECULATE instance Additive.C (T Double) #-}     {-# INLINE zero #-}-    zero                        =  Cons zero zero+    {-# INLINE negate #-}     {-# INLINE (+) #-}-    (Cons x y) + (Cons x' y')   =  Cons (x+x') (y+y')     {-# INLINE (-) #-}-    (Cons x y) - (Cons x' y')   =  Cons (x-x') (y-y')-    {-# INLINE negate #-}-    negate (Cons x y)           =  Cons (negate x) (negate y)+    zero   = Cons zero zero+    (+)    = Elem.run2 $ Elem.with Cons <*>.+  real <*>.+  imag+    (-)    = Elem.run2 $ Elem.with Cons <*>.-  real <*>.-  imag+    negate = Elem.run  $ Elem.with Cons <*>.-$ real <*>.-$ imag  instance  (Ring.C a) => Ring.C (T a)  where     {-# SPECULATE instance Ring.C (T Float) #-}@@ -295,7 +319,8 @@ --   because it requires the Algebra.Module constraint instance (Module.C a b) => Module.C a (T b) where    {-# INLINE (*>) #-}-   s *> (Cons x y)  = Cons (s *> x) (s *> y)+   (*>) = Elem.run2 $ Elem.with Cons <*>.*> real <*>.*> imag+   -- s *> (Cons x y)  = Cons (s *> x) (s *> y)  instance (VectorSpace.C a b) => VectorSpace.C a (T b)
src/Number/Quaternion.hs view
@@ -47,13 +47,15 @@ import qualified Algebra.ZeroTestable as ZeroTestable  import Algebra.ZeroTestable(isZero)-import Algebra.Module((*>))--- import Algebra.Algebraic((^/))+import Algebra.Module((*>), (<*>.*>), )  import qualified Number.Complex as Complex  import Number.Complex ((+:)) +import qualified NumericPrelude.Elementwise as Elem+import Algebra.Additive ((<*>.+), (<*>.-), (<*>.-$), )+ -- import qualified Data.Typeable as Ty import Data.Array (Array, (!)) import qualified Data.Array as Array@@ -256,10 +258,10 @@ instance (Additive.C a) => Additive.C (T a)  where    {-# SPECIALISE instance Additive.C (T Float) #-}    {-# SPECIALISE instance Additive.C (T Double) #-}-   zero			=  Cons zero zero-   (Cons xr xi) + (Cons yr yi)	=  Cons (xr+yr) (xi+yi)-   (Cons xr xi) - (Cons yr yi)	=  Cons (xr-yr) (xi-yi)-   negate (Cons x y)		=  Cons (negate x) (negate y)+   zero   = Cons zero zero+   (+)    = Elem.run2 $ Elem.with Cons <*>.+  real <*>.+  imag+   (-)    = Elem.run2 $ Elem.with Cons <*>.-  real <*>.-  imag+   negate = Elem.run  $ Elem.with Cons <*>.-$ real <*>.-$ imag  instance (Ring.C a) => Ring.C (T a)  where    {-# SPECIALISE instance Ring.C (T Float) #-}@@ -288,7 +290,7 @@ -- | The '(*>)' method can't replace 'scale' --   because it requires the Algebra.Module constraint instance (Module.C a b) => Module.C a (T b) where-   s *> (Cons r i)  = Cons (s *> r) (s *> i)+   (*>) = Elem.run2 $ Elem.with Cons <*>.*> real <*>.*> imag  instance (VectorSpace.C a b) => VectorSpace.C a (T b)
src/Number/Ratio.hs view
@@ -40,6 +40,10 @@  import Control.Monad(liftM, liftM2, ) +import Foreign.Storable (Storable (..), )+import qualified Foreign.Storable.Record as Store+import Control.Applicative (liftA2, )+ import Test.QuickCheck (Arbitrary(arbitrary,coarbitrary)) import System.Random (Random(..), RandomGen, ) @@ -157,6 +161,22 @@       liftM2 (%) arbitrary          (liftM (\x -> if isZero x then one else x) arbitrary)    coarbitrary = undefined+++instance (Storable a, PID.C a) => Storable (T a) where+   sizeOf    = Store.sizeOf store+   alignment = Store.alignment store+   peek      = Store.peek store+   poke      = Store.poke store++store ::+   (Storable a, PID.C a) =>+   Store.Dictionary (T a)+store =+   Store.run $+   liftA2 (%)+      (Store.element numerator)+      (Store.element denominator)  {- This instance may not be appropriate for mathematical objects other than numbers.
+ src/NumericPrelude/Elementwise.hs view
@@ -0,0 +1,41 @@+module NumericPrelude.Elementwise where++import Control.Applicative (Applicative(pure, (<*>)), )++{- |+A reader monad for the special purpose+of defining instances of certain operations on tuples and records.+It does not add any new functionality to the common Reader monad,+but it restricts the functions to the required ones+and exports them from one module.+That is you do not have to import+both Control.Monad.Trans.Reader and Control.Applicative.+The type also tells the user, for what the Reader monad is used.+We can more easily replace or extend the implementation when needed.+-}+newtype T v a = Cons {run :: v -> a}++{-# INLINE with #-}+with :: a -> T v a+with e = Cons $ \ _v -> e++{-# INLINE element #-}+element :: (v -> a) -> T v a+element = Cons+++{-# INLINE run2 #-}+run2 :: T (x,y) a -> x -> y -> a+run2 = curry . run++instance Functor (T v) where+   {-# INLINE fmap #-}+   fmap f (Cons e) =+      Cons $ \v -> f $ e v++instance Applicative (T v) where+   {-# INLINE pure #-}+   {-# INLINE (<*>) #-}+   pure = with+   Cons f <*> Cons e =+      Cons $ \v -> f v $ e v
+ test/Gaussian.hs view
@@ -0,0 +1,6 @@+module Main where++import qualified MathObj.Gaussian.Example as Example++main :: IO ()+main = Example.polyApprox
test/Test/MathObj/Gaussian/Bell.hs view
@@ -62,6 +62,8 @@           simple $ \x -> nest 2 G.fourier x == G.reverse x) :       ("reverse identity",           simple $ \x -> nest 2 G.reverse x == x) :+      ("fourier unit",+          quickCheck $ G.fourier G.unit == (G.unit :: G.T Rational)) :       ("translate additive",           simple $ \x a b ->              G.translate a (G.translate b x) == G.translate (a+b) x) :@@ -90,5 +92,5 @@              G.shrink a x == G.dilate (recip a) x) :       ("fourier dilate",           simple $ \x a -> a>0 ==>-             G.fourier (G.dilate a x) == G.shrink a (G.fourier x)) :+             G.fourier (G.dilate a x) == G.amplify a (G.shrink a (G.fourier x))) :       []
test/Test/MathObj/Gaussian/Variance.hs view
@@ -58,5 +58,5 @@              G.shrink a x == G.dilate (recip a) x) :       ("fourier dilate",           simple $ \x a -> a>0 ==>-             G.fourier (G.dilate a x) == G.shrink a (G.fourier x)) :+             G.fourier (G.dilate a x) == G.amplify a (G.shrink a (G.fourier x))) :       []
test/Test/Run.hs view
@@ -1,5 +1,6 @@ module Main where +import qualified Test.MathObj.Gaussian.Polynomial as GaussPoly import qualified Test.MathObj.Gaussian.Variance as GaussVariance import qualified Test.MathObj.Gaussian.Bell as GaussBell import qualified Test.MathObj.PartialFraction as PartialFraction@@ -13,6 +14,7 @@    do HUnitText.runTestTT (HUnit.TestList $          GaussVariance.tests :          GaussBell.tests :+         GaussPoly.tests :          PartialFraction.tests :          Polynomial.tests :          PowerSeries.tests :