packages feed

interpolation 0.1.0.3 → 0.1.1

raw patch · 4 files changed

+164/−82 lines, 4 filesdep +comfort-arraydep +lapackdep +tfpdep −hmatrixdep −hmatrix-bandeddep ~gnuplotPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: comfort-array, lapack, tfp

Dependencies removed: hmatrix, hmatrix-banded

Dependency ranges changed: gnuplot

API changes (from Hackage documentation)

- Numeric.Interpolation.Basis: coefficientsToCubicLinear :: (Fractional a) => [a] -> [a] -> T a (a, a)
+ Numeric.Interpolation.Basis: coefficientsToCubicLinear :: Fractional a => [a] -> [a] -> T a (a, a)
- Numeric.Interpolation.Basis: coefficientsToCubicParabola :: (Fractional a) => [a] -> [a] -> T a (a, a)
+ Numeric.Interpolation.Basis: coefficientsToCubicParabola :: Fractional a => [a] -> [a] -> T a (a, a)
- Numeric.Interpolation.Basis: cubicLinear :: (Fractional a) => [a] -> [T a (a, a)]
+ Numeric.Interpolation.Basis: cubicLinear :: Fractional a => [a] -> [T a (a, a)]
- Numeric.Interpolation.Basis: cubicParabola :: (Fractional a) => [a] -> [T a (a, a)]
+ Numeric.Interpolation.Basis: cubicParabola :: Fractional a => [a] -> [T a (a, a)]
- Numeric.Interpolation.Basis: hermite1 :: (Num b) => [a] -> [T a (b, b)]
+ Numeric.Interpolation.Basis: hermite1 :: Num b => [a] -> [T a (b, b)]
- Numeric.Interpolation.Basis: linear :: (Num b) => [a] -> [T a b]
+ Numeric.Interpolation.Basis: linear :: Num b => [a] -> [T a b]
- Numeric.Interpolation.Basis.Compact: cubicLinear :: (Fractional a) => [a] -> [T a (a, a)]
+ Numeric.Interpolation.Basis.Compact: cubicLinear :: Fractional a => [a] -> [T a (a, a)]
- Numeric.Interpolation.Basis.Compact: cubicParabola :: (Fractional a) => [a] -> [T a (a, a)]
+ Numeric.Interpolation.Basis.Compact: cubicParabola :: Fractional a => [a] -> [T a (a, a)]
- Numeric.Interpolation.Basis.Compact: hermite1 :: (Num b) => [a] -> [T a (b, b)]
+ Numeric.Interpolation.Basis.Compact: hermite1 :: Num b => [a] -> [T a (b, b)]
- Numeric.Interpolation.Basis.Compact: linear :: (Num b) => [a] -> [T a b]
+ Numeric.Interpolation.Basis.Compact: linear :: Num b => [a] -> [T a b]
- Numeric.Interpolation.Basis.Full: hermite1 :: (Num b) => [a] -> [T a (b, b)]
+ Numeric.Interpolation.Basis.Full: hermite1 :: Num b => [a] -> [T a (b, b)]
- Numeric.Interpolation.Basis.Full: linear :: (Num b) => [a] -> [T a b]
+ Numeric.Interpolation.Basis.Full: linear :: Num b => [a] -> [T a b]
- Numeric.Interpolation.NodeList: Node :: (x, y) -> (T x y) -> (T x y) -> T x y
+ Numeric.Interpolation.NodeList: Node :: (x, y) -> T x y -> T x y -> T x y
- Numeric.Interpolation.Piece: hermite1 :: (Fractional a) => T a a (a, a)
+ Numeric.Interpolation.Piece: hermite1 :: Fractional a => T a a (a, a)
- Numeric.Interpolation.Piece: linear :: (Fractional a) => T a a a
+ Numeric.Interpolation.Piece: linear :: Fractional a => T a a a
- Numeric.Interpolation.Piecewise: interpolate :: (Ord x) => T x y ny -> T x ny -> x -> y
+ Numeric.Interpolation.Piecewise: interpolate :: Ord x => T x y ny -> T x ny -> x -> y
- Numeric.Interpolation.Piecewise: interpolateConstantExt :: (Ord x) => T x y ny -> T x ny -> x -> y
+ Numeric.Interpolation.Piecewise: interpolateConstantExt :: Ord x => T x y ny -> T x ny -> x -> y
- Numeric.Interpolation.Type: Cons :: [x] -> [y] -> String -> T x y ny -> Int -> [x] -> [T x ny] -> [x] -> x -> [(Int, y)] -> [x] -> [y] -> T x ny -> ny -> y -> T x y ny
+ Numeric.Interpolation.Type: Cons :: ([x] -> [y] -> String) -> T x y ny -> Int -> ([x] -> [T x ny]) -> ([x] -> x -> [(Int, y)]) -> ([x] -> [y] -> T x ny) -> (ny -> y) -> T x y ny

Files

example/Fit.hs view
@@ -1,31 +1,50 @@+{-# LANGUAGE Rank2Types #-} module Main where +import qualified Plot2DExtra+ import qualified Numeric.Interpolation.NodeList as Nodes import qualified Numeric.Interpolation.Piecewise as Piecewise import qualified Numeric.Interpolation.Type as Type -import qualified Data.Packed.ST as PackST-import qualified Data.Packed.Matrix as Matrix-import qualified Data.Packed.Vector as Vector-import Data.Packed.Matrix (Matrix)-import Data.Packed.Vector (Vector)+import qualified Numeric.LAPACK.Singular as Singular+import qualified Numeric.LAPACK.ShapeStatic as ShapeStatic+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape+import qualified Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite+                                                       as BandedSPD+import qualified Numeric.LAPACK.Matrix.BandedHermitian as BandedHermitian+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector -import qualified Numeric.LinearAlgebra.Banded as Banded-import qualified Numeric.Container as Container-import Numeric.Container ((<\>))+import qualified Type.Data.Num.Unary as Unary+import Type.Base.Proxy (Proxy(Proxy)) +import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+ import qualified Graphics.Gnuplot.Advanced as GP+import qualified Graphics.Gnuplot.Frame as Frame+import qualified Graphics.Gnuplot.Frame.OptionSet as Opts import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D+import qualified Graphics.Gnuplot.LineSpecification as LineSpec  import System.Random (randomRs, mkStdGen) import Control.Monad.HT (void)-import Control.Monad (when, zipWithM_, forM_)+import Control.Applicative ((<$>))  import qualified Data.Foldable as Fold+import Data.Tuple.HT (mapFst) import Data.Monoid ((<>))  +type BandedHermitianMatrix k =+       BandedHermitian.BandedHermitian k Matrix.ZeroInt+type Matrix = Matrix.General Matrix.ZeroInt Matrix.ZeroInt+type Vector = Vector.Vector Matrix.ZeroInt+type ShortVector k = Vector.Vector (ShapeStatic.ZeroBased k)++ noisy :: [(Double, Double)] noisy =    take 100 $@@ -37,19 +56,24 @@ basisMatrixFull ::    Type.T Double Double ny -> [Double] -> [Double] -> Matrix Double basisMatrixFull typ xs txs0 =-   let txs = Vector.fromList txs0-   in  Matrix.fromColumns $-       map (flip Container.cmap txs . Piecewise.interpolateConstantExt typ) $+   let txs = Vector.autoFromList txs0+   in  Matrix.fromColumns (Array.shape txs) $+       map (flip Array.map txs . Piecewise.interpolateConstantExt typ) $        Type.basisFunctions typ xs +zipRowsWith :: (a -> b -> c) -> [a] -> [[b]] -> [c]+zipRowsWith f as bs = concat $ zipWith (map . f) as bs++-- ToDo: generalize shapes basisMatrixSparse ::    Type.T Double Double ny -> [Double] -> [Double] -> Matrix Double-basisMatrixSparse typ xs txs = PackST.runSTMatrix $ do-   mat <- PackST.newMatrix 0 (length txs) (length $ Type.basisFunctions typ xs)-   zipWithM_-      (\k -> mapM_ (uncurry (PackST.writeMatrix mat k))) [0..] $-      map (Type.sampleBasisFunctions typ xs) txs-   return mat+basisMatrixSparse typ xs txs =+   Array.fromAssociations+      (MatrixShape.general MatrixShape.RowMajor+         (Shape.ZeroBased $ length txs)+         (Shape.ZeroBased $ length $ Type.basisFunctions typ xs)) 0 $+   zipRowsWith (\k (j,x) -> ((k,j),x)) [0..] $+   map (Type.sampleBasisFunctions typ xs) txs  fit ::    Type.T Double Double ny ->@@ -58,39 +82,56 @@    let (txs, tys) = unzip target        matrix = basisMatrixSparse typ xs txs    in  Type.coefficientsToInterpolator typ xs $-       Vector.toList $ matrix <\> Vector.fromList tys+       Vector.toList $+       Matrix.unliftColumn MatrixShape.ColumnMajor+         (snd . Singular.leastSquaresMinimumNormRCond 1e-5 matrix) $+       Vector.autoFromList tys  matrixDiff ::-   Type.T Double Double ny ->-   [Double] -> [(Double, Double)] -> Double+   Type.T Double Double ny -> [Double] -> [(Double, Double)] -> Double matrixDiff typ xs target =    let txs = map fst target-   in  Container.maxElement $ Container.cmap abs $-       Container.sub+   in  snd $ Vector.argAbs1Maximum $+       Matrix.sub           (basisMatrixFull typ xs txs)           (basisMatrixSparse typ xs txs)   mulSparseMatrixVector ::    Int -> [[(Int, Double)]] -> [Double] -> Vector Double-mulSparseMatrixVector size samples tys = PackST.runSTVector $ do-   vec <- PackST.newVector 0 size-   forM_ (zip samples tys) $ \(row,ty) ->-      forM_ row $ \(k,y) ->-         PackST.modifyVector vec k (+y*ty)-   return vec+mulSparseMatrixVector size samples tys =+   Array.accumulate (+)+      (Vector.constant (Matrix.zeroInt size) 0)+      (zipRowsWith (\ty (k,y) -> (k, y*ty)) tys samples) +{- ToDo:+use index type (ZeroInt, Enumeration Order)+with Order = Absolute | Derivative+Problem: not all interpolation types use derivatives+-}+shortVector ::+   (Unary.Natural k) =>+   Proxy k -> [(Int, Double)] -> (Int, ShortVector k Double)+shortVector width xs =+   let i0 = minimum $ map fst xs+   in (i0,+--       Array.fromAssociations (ShapeStatic.ZeroBased Proxy) 0 $+       Array.reshape (ShapeStatic.ZeroBased Proxy) $+       Array.fromAssociations (Matrix.zeroInt $ Unary.integralFromProxy width) 0 $+       map (mapFst (subtract i0)) xs)+ bandedGramian ::-   Int -> Int -> [[(Int, Double)]] -> Banded.SymmetricMatrix Double+   (Unary.Natural k) =>+   Int -> Proxy (Unary.Succ k) ->+   [[(Int, Double)]] -> BandedHermitianMatrix k Double bandedGramian size width samples =-      Banded.SymmetricMatrix $ PackST.runSTMatrix $ do-   mat <- PackST.newMatrix 0 size width-   forM_ samples $ \row ->-      forM_ row $ \(k,yk) ->-      forM_ row $ \(j,yj) ->-         when (k<=j) $ PackST.modifyMatrix mat k (j-k) (+yk*yj)-   return mat+   BandedHermitian.sumRank1 MatrixShape.ColumnMajor (Matrix.zeroInt size) $+   map ((,) 1 . shortVector width) samples +reifyPositive ::+   Integer -> (forall s. Unary.Natural s => Proxy (Unary.Succ s) -> w) -> w+reifyPositive n f = Unary.reifyNatural (n-1) (f . Unary.succ)+ fitBanded ::    Type.T Double Double ny ->    [Double] -> [(Double, Double)] -> Nodes.T Double ny@@ -98,11 +139,12 @@    let size = length $ Type.basisFunctions typ xs        (txs, tys) = unzip target        samples = map (Type.sampleBasisFunctions typ xs) txs-       matrix =-          Banded.choleskyDecompose $-          bandedGramian size (Type.basisOverlap typ) samples-   in  Type.coefficientsToInterpolator typ xs $ Vector.toList $-       Banded.choleskySolve matrix $ mulSparseMatrixVector size samples tys+   in reifyPositive (toInteger $ Type.basisOverlap typ) +         (\width ->+            Type.coefficientsToInterpolator typ xs $ Vector.toList $+            Matrix.unliftColumn MatrixShape.ColumnMajor+               (BandedSPD.solve (bandedGramian size width samples)) $+            mulSparseMatrixVector size samples tys)  bandedDiff ::    (ny -> ny -> Double) ->@@ -122,44 +164,62 @@   plotBasisFunctions ::-   Type.T Double Double ny -> [Double] -> Plot2D.T Double Double-plotBasisFunctions nodeType xs =+   (String, Type.T Double Double ny) -> [Double] ->+   Frame.T (Graph2D.T Double Double)+plotBasisFunctions (title, nodeType) xs =    let abscissa = Plot2D.linearScale 1000 (minimum xs, maximum xs)-   in  Plot2D.functions Graph2D.lines abscissa $+   in  Frame.cons (Opts.title title $ Opts.key False Opts.deflt) $+       Plot2D.functions Graph2D.lines abscissa $        map (Piecewise.interpolateConstantExt nodeType) $        Type.basisFunctions nodeType xs  +typeLinear :: (String, Type.T Double Double Double)+typeLinear        = ("linear",        Type.linear)+typeHermite1, typeCubicLinear, typeCubicParabola ::+   (String, Type.T Double Double (Double, Double))+typeHermite1      = ("hermite1",      Type.hermite1)+typeCubicLinear   = ("cubicLinear",   Type.cubicLinear)+typeCubicParabola = ("cubicParabola", Type.cubicParabola)++ main :: IO () main = do    let xs = [0, 1, 3, 4, 6, 7]        exs = (-1) : xs ++ [8]-   void $ GP.plotDefault $ plotBasisFunctions Type.linear xs-   void $ GP.plotDefault $ plotBasisFunctions Type.hermite1 xs-   void $ GP.plotDefault $ plotBasisFunctions Type.cubicLinear exs-   void $ GP.plotDefault $ plotBasisFunctions Type.cubicParabola exs+   void $ GP.plotDefault $ plotBasisFunctions typeLinear xs+   void $ GP.plotDefault $ plotBasisFunctions typeHermite1 xs+   void $ GP.plotDefault $ plotBasisFunctions typeCubicLinear exs+   void $ GP.plotDefault $ plotBasisFunctions typeCubicParabola exs    let linearNodes = fit Type.linear xs noisy        hermite1Nodes = fit Type.hermite1 xs noisy        cubicLinearNodes = fit Type.cubicLinear exs noisy        cubicParabolaNodes = fit Type.cubicParabola exs noisy    void $ GP.plotDefault $-      Plot2D.list Graph2D.points noisy+      (Graph2D.lineSpec (LineSpec.title "target" LineSpec.deflt)+         <$> Plot2D.list Graph2D.points noisy)       <>-      (Plot2D.functions Graph2D.lines (Plot2D.linearScale 1000 (-2,10)) $-       Piecewise.interpolateConstantExt Type.linear linearNodes :-       Piecewise.interpolateConstantExt Type.hermite1 hermite1Nodes :-       Piecewise.interpolateConstantExt Type.cubicLinear cubicLinearNodes :-       Piecewise.interpolateConstantExt Type.cubicParabola cubicParabolaNodes :-       [])+      (let interpolate (name,typ) nodes =+             (name, Piecewise.interpolateConstantExt typ nodes)+       in Plot2DExtra.functions (Plot2D.linearScale 1000 (-2,10)) $+            interpolate typeLinear linearNodes :+            interpolate typeHermite1 hermite1Nodes :+            interpolate typeCubicLinear cubicLinearNodes :+            interpolate typeCubicParabola cubicParabolaNodes :+            [])     putStrLn "differences between matrices should be almost zero:"-   putStrLn $ "linear: " ++ show (matrixDiff Type.linear xs noisy)-   putStrLn $ "hermite1: " ++ show (matrixDiff Type.hermite1 xs noisy)-   putStrLn $ "cubicLinear: " ++ show (matrixDiff Type.cubicLinear exs noisy)-   putStrLn $ "cubicParabola: " ++ show (matrixDiff Type.cubicParabola exs noisy)+   let printMatrixDiff (name,typ) ps =+         putStrLn $ name ++ ": " ++ show (matrixDiff typ ps noisy)+   printMatrixDiff typeLinear xs+   printMatrixDiff typeHermite1 xs+   printMatrixDiff typeCubicLinear exs+   printMatrixDiff typeCubicParabola exs     putStrLn "differences between samples should be almost zero:"-   putStrLn $ "linear: " ++ show (bandedDiff absDiffSingle Type.linear xs noisy)-   putStrLn $ "hermite1: " ++ show (bandedDiff absDiffPair Type.hermite1 xs noisy)-   putStrLn $ "cubicLinear: " ++ show (bandedDiff absDiffPair Type.cubicLinear exs noisy)-   putStrLn $ "cubicParabola: " ++ show (bandedDiff absDiffPair Type.cubicParabola exs noisy)+   let printBandedDiff diff (name,typ) ps =+         putStrLn $ name ++ ": " ++ show (bandedDiff diff typ ps noisy)+   printBandedDiff absDiffSingle typeLinear xs+   printBandedDiff absDiffPair typeHermite1 xs+   printBandedDiff absDiffPair typeCubicLinear exs+   printBandedDiff absDiffPair typeCubicParabola exs
example/Plot.hs view
@@ -1,5 +1,7 @@ module Main where +import qualified Plot2DExtra+ import qualified Numeric.Interpolation.NodeList as Nodes import qualified Numeric.Interpolation.Piecewise as Piecewise import qualified Numeric.Interpolation.Basis as Basis@@ -9,7 +11,6 @@ import qualified Graphics.Gnuplot.Advanced as GP  import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D-import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D  import Control.Monad.HT (void) @@ -24,10 +25,11 @@        cubicLinearNodes = Basis.coefficientsToCubicLinear xs $ map sin xs        cubicParabolaNodes = Basis.coefficientsToCubicParabola xs $ map sin xs    void $ GP.plotDefault $-      Plot2D.functions Graph2D.lines-         (Plot2D.linearScale 1000 (-2,15))-         [Piecewise.interpolateConstantExt Type.linear linearNodes,-          Piecewise.interpolateConstantExt Type.hermite1 hermite1Nodes,-          Piecewise.interpolateConstantExt Type.hermite1 cubicLinearNodes,-          Piecewise.interpolateConstantExt Type.hermite1 cubicParabolaNodes,-          sin]+      Plot2DExtra.functions (Plot2D.linearScale 1000 (-2,15)) $+         let interpolate typ nodes = Piecewise.interpolateConstantExt typ nodes+         in ("linear",        interpolate Type.linear linearNodes) :+            ("hermite1",      interpolate Type.hermite1 hermite1Nodes) :+            ("cubicLinear",   interpolate Type.hermite1 cubicLinearNodes) :+            ("cubicParabola", interpolate Type.hermite1 cubicParabolaNodes) :+            ("sin", sin) :+            []
+ example/Plot2DExtra.hs view
@@ -0,0 +1,17 @@+module Plot2DExtra where++import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D+import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D+import qualified Graphics.Gnuplot.LineSpecification as LineSpec+import qualified Graphics.Gnuplot.Value.Tuple as Tuple+import qualified Graphics.Gnuplot.Value.Atom as Atom++import Data.Tuple.HT (mapFst)+++functions ::+   (Atom.C x, Atom.C y, Tuple.C x, Tuple.C y) =>+   [x] -> [(String, x -> y)] -> Plot2D.T x y+functions scale =+   Plot2D.functionsWithLineSpec Graph2D.lines scale .+   map (mapFst (flip LineSpec.title LineSpec.deflt))
interpolation.cabal view
@@ -1,5 +1,5 @@ Name:             interpolation-Version:          0.1.0.3+Version:          0.1.1 License:          BSD3 License-File:     LICENSE Author:           Henning Thielemann@@ -26,7 +26,7 @@   * @example/Fit.hs@:     Demonstrates how to use the basis functions     for fitting an interpolation function to a given function-    using a linear least squares solver like @<\>@ from @hmatrix@.+    using a linear least squares solver like @<\>@ from @lapack@.     We use a distorted sinus as target.   .   The package needs only Haskell 98.@@ -42,12 +42,12 @@   description: Build example executables   default:     False -Flag hmatrix-  description: Build examples that depend on hmatrix+Flag lapack+  description: Build examples that depend on lapack   default:     True  Source-Repository this-  Tag:         0.1.0.3+  Tag:         0.1.1   Type:        darcs   Location:    http://hub.darcs.net/thielema/interpolation/ @@ -78,12 +78,13 @@  Executable interpolation-plot   Main-Is:             Plot.hs+  Other-Modules:       Plot2DExtra   Hs-Source-Dirs:      example   GHC-Options:         -Wall   If flag(buildExamples)     Build-Depends:       interpolation,-      gnuplot >=0.5.2 && <0.6,+      gnuplot >=0.5.6 && <0.6,       utility-ht >=0.0.9 && <0.1,       base >=4.5 && <5   Else@@ -91,15 +92,17 @@  Executable interpolation-fit   Main-Is:             Fit.hs+  Other-Modules:       Plot2DExtra   Hs-Source-Dirs:      example   GHC-Options:         -Wall-  If flag(buildExamples) && flag(hmatrix)+  If flag(buildExamples) && flag(lapack)     Build-Depends:       interpolation,-      hmatrix-banded >=0.0 && <0.1,-      hmatrix >=0.15 && <0.16,+      lapack >=0.2.2 && <0.3,+      tfp >=1.0 && <1.1,+      comfort-array >=0.3 && <0.4,       random >=1.0 && <1.2,-      gnuplot >=0.5.2 && <0.6,+      gnuplot >=0.5.6 && <0.6,       utility-ht >=0.0.9 && <0.1,       base >=4.5 && <5   Else