packages feed

interpolation (empty) → 0.0

raw patch · 16 files changed

+773/−0 lines, 16 filesdep +QuickCheckdep +basedep +gnuplotsetup-changed

Dependencies added: QuickCheck, base, gnuplot, hmatrix, interpolation, random, utility-ht

Files

+ LICENSE view
@@ -0,0 +1,26 @@+Copyright (c) Henning Thielemann 2014+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:+1. Redistributions of source code must retain the above copyright+   notice, this list of conditions and the following disclaimer.+2. Redistributions in binary form must reproduce the above copyright+   notice, this list of conditions and the following disclaimer in the+   documentation and/or other materials provided with the distribution.+3. Neither the name of the University nor the names of its contributors+   may be used to endorse or promote products derived from this software+   without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+SUCH DAMAGE.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ example/Fit.hs view
@@ -0,0 +1,71 @@+module Main where++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.Matrix as Matrix+import qualified Data.Packed.Vector as Vector+import qualified Numeric.Container as Container+import Numeric.Container ((<\>))++import qualified Graphics.Gnuplot.Advanced as GP+import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D+import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D++import System.Random (randomRs, mkStdGen)+import Control.Monad.HT (void)+import Data.Monoid ((<>))+++noisy :: [(Double, Double)]+noisy =+   take 100 $+   zipWith+      (\x d -> (x, sin x + d))+      (randomRs (0,7) (mkStdGen 23))+      (randomRs (-0.2,0.2) (mkStdGen 42))++fit ::+   Type.T Double Double ny ->+   [Double] -> [(Double, Double)] -> Nodes.T Double ny+fit typ xs target =+   let txs = Vector.fromList $ map fst target+       tys = Vector.fromList $ map snd target+       matrix =+          Matrix.fromColumns $+          map (flip Container.cmap txs . Piecewise.interpolateConstantExt typ) $+          Type.basisFunctions typ xs+   in  Type.coefficientsToInterpolator typ xs $+       Vector.toList $ matrix <\> tys++plotBasisFunctions ::+   Type.T Double Double ny -> [Double] -> Plot2D.T Double Double+plotBasisFunctions nodeType xs =+   let abscissa = Plot2D.linearScale 1000 (minimum xs, maximum xs)+   in  Plot2D.functions Graph2D.lines abscissa $+       map (Piecewise.interpolateConstantExt nodeType) $+       Type.basisFunctions nodeType xs+++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.cubic xs+   void $ GP.plotDefault $ plotBasisFunctions Type.cubicLinear exs+   void $ GP.plotDefault $ plotBasisFunctions Type.cubicParabola exs+   let linearNodes = fit Type.linear xs noisy+       hermite1Nodes = fit Type.cubic xs noisy+       cubicLinearNodes = fit Type.cubicLinear exs noisy+       cubicParabolaNodes = fit Type.cubicParabola exs noisy+   void $ GP.plotDefault $+      Plot2D.list Graph2D.points noisy+      <>+      (Plot2D.functions Graph2D.lines (Plot2D.linearScale 1000 (-2,10)) $+       Piecewise.interpolateConstantExt Type.linear linearNodes :+       Piecewise.interpolateConstantExt Type.cubic hermite1Nodes :+       Piecewise.interpolateConstantExt Type.cubicLinear cubicLinearNodes :+       Piecewise.interpolateConstantExt Type.cubicParabola cubicParabolaNodes :+       [])
+ example/Plot.hs view
@@ -0,0 +1,33 @@+module Main where++import qualified Numeric.Interpolation.NodeList as Nodes+import qualified Numeric.Interpolation.Piecewise as Piecewise+import qualified Numeric.Interpolation.Basis as Basis+import qualified Numeric.Interpolation.Type as Type+++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)+++xs :: [Double]+xs = [0, 1, 3, 4, 6, 7, 9, 10, 11, 13]++main :: IO ()+main = do+   let linearNodes = Nodes.fromList $ map (\x -> (x, sin x)) xs+       hermite1Nodes = Nodes.fromList $ map (\x -> (x, (sin x, cos x))) xs+       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.cubic hermite1Nodes,+          Piecewise.interpolateConstantExt Type.cubic cubicLinearNodes,+          Piecewise.interpolateConstantExt Type.cubic cubicParabolaNodes,+          sin]
+ interpolation.cabal view
@@ -0,0 +1,112 @@+Name:             interpolation+Version:          0.0+License:          BSD3+License-File:     LICENSE+Author:           Henning Thielemann+Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>+Homepage:         http://code.haskell.org/~thielema/interpolation/+Category:         Math+Synopsis:         piecewise linear and cubic Hermite interpolation+Description:+  Represent real functions by linear or cubic segments.+  The package provides both data structures+  for efficient lookup of interpolation intervals,+  and computation of basis functions.+  .+  There are two examples that can be built with+  .+  > cabal install -fbuildExamples+  .+  * @example/Plot.hs@:+    Interpolate a sinus curve using piecewise linear interpolation+    and piecewise Hermite cubic interpolation.+    For the latter one we provide the derivatives of the sinus function+    at the interpolation nodes.+  .+  * @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@.+    We use a distorted sinus as target.+  .+  The package needs only Haskell 98.+Tested-With:      GHC==7.4.2, GHC==7.6.3, GHC==7.8.2+Cabal-Version:    >=1.8+Build-Type:       Simple++Flag buildExamples+  description: Build example executables+  default:     False++Flag hmatrix+  description: Build examples that depend on hmatrix+  default:     True++Source-Repository this+  Tag:         0.0+  Type:        darcs+  Location:    http://code.haskell.org/~thielema/interpolation/++Source-Repository head+  Type:        darcs+  Location:    http://code.haskell.org/~thielema/interpolation/++Library+  Build-Depends:+    utility-ht >=0.0.1 && <0.1,+    base >=4 && <5++  GHC-Options:      -Wall+  Hs-Source-Dirs:   src, private+  Exposed-Modules:+    Numeric.Interpolation.NodeList+    Numeric.Interpolation.Piece+    Numeric.Interpolation.Piecewise+    Numeric.Interpolation.Type+    Numeric.Interpolation.Basis+    Numeric.Interpolation.Basis.Compact+    Numeric.Interpolation.Basis.Full+  Other-Modules:+    Numeric.Interpolation.Private.Piece+    Numeric.Interpolation.Private.Basis++Executable interpolation-plot+  Main-Is:             Plot.hs+  Hs-Source-Dirs:      example+  GHC-Options:         -Wall+  If flag(buildExamples)+    Build-Depends:+      interpolation,+      gnuplot >=0.5.2 && <0.6,+      utility-ht >=0.0.9 && <0.1,+      base >=4.5 && <4.8+  Else+    Buildable: False++Executable interpolation-fit+  Main-Is:             Fit.hs+  Hs-Source-Dirs:      example+  GHC-Options:         -Wall+  If flag(buildExamples) && flag(hmatrix)+    Build-Depends:+      interpolation,+      hmatrix >=0.15 && <0.16,+      random >=1.0 && <1.1,+      gnuplot >=0.5.2 && <0.6,+      utility-ht >=0.0.9 && <0.1,+      base >=4.5 && <4.8+  Else+    Buildable: False++Test-Suite interpolation-test+  Type:                exitcode-stdio-1.0+  Main-Is:             Test.hs+  Other-Modules:+    Test.Piece+  Hs-Source-Dirs:      test, private+  GHC-Options:         -Wall+  Build-Depends:+    interpolation,+    QuickCheck >=2.4 && <2.8,+    utility-ht >=0.0.9 && <0.1,+    base >=4.5 && <4.8
+ private/Numeric/Interpolation/Private/Basis.hs view
@@ -0,0 +1,53 @@+module Numeric.Interpolation.Private.Basis where++import Numeric.Interpolation.Private.Piece (sqr)++import qualified Data.List.Match as Match+++hermite1Split :: [a] -> [b] -> [(b, b)]+hermite1Split xs = uncurry zip . Match.splitAt xs+++parabolaDerivative ::+   (Fractional a) => (a,a) -> (a,a) -> (a,a) -> a -> (a,a)+parabolaDerivative (x0,y0) (x1,y1) (x2,y2) x =+   let l0 = (x-x1)*(x-x2)/((x0-x1)*(x0-x2))+       l1 = (x-x0)*(x-x2)/((x1-x0)*(x1-x2))+       l2 = (x-x0)*(x-x1)/((x2-x0)*(x2-x1))+       dl0 = (2*x-x1-x2)/((x0-x1)*(x0-x2))+       dl1 = (2*x-x0-x2)/((x1-x0)*(x1-x2))+       dl2 = (2*x-x0-x1)/((x2-x0)*(x2-x1))+   in  (y0*l0 + y1*l1 + y2*l2, y0*dl0 + y1*dl1 + y2*dl2)++parabolaBasisDerivativeLeft,+   parabolaBasisDerivativeCenter,+   parabolaBasisDerivativeRight ::+      (Fractional a) => a -> a -> a -> a+parabolaBasisDerivativeLeft   x0 x1 x2 = (x1-x2)/((x0-x1)*(x0-x2))+parabolaBasisDerivativeCenter x0 x1 x2 = 1/(x1-x0) + 1/(x1-x2)+parabolaBasisDerivativeRight  x0 x1 x2 = (x1-x0)/((x2-x0)*(x2-x1))++parabolaDerivativeCenterNode ::+   (Fractional a) => (a,a) -> (a,a) -> (a,a) -> a+parabolaDerivativeCenterNode (x0,y0) (x1,y1) (x2,y2) =+   y0 * parabolaBasisDerivativeLeft   x0 x1 x2 ++   y1 * parabolaBasisDerivativeCenter x0 x1 x2 ++   y2 * parabolaBasisDerivativeRight  x0 x1 x2+++parabola2ndDerivativeCenterNode ::+   (Fractional a) => (a,a) -> (a,a) -> (a,a) -> (a,a) -> a+parabola2ndDerivativeCenterNode (xl,yl) (x0,y0) (x1,y1) (x2,y2) =+   let dy0 =+          yl * (x0-x1)/((xl-x0)*(xl-x1)) ++          y0 * (1/(x0-xl) + 1/(x0-x1)) ++          y1 * (x0-xl)/((x1-xl)*(x1-x0))+       dy1 =+          y0 * (x1-x2)/((x0-x1)*(x0-x2)) ++          y1 * (1/(x1-x0) + 1/(x1-x2)) ++          y2 * (x1-x0)/((x2-x0)*(x2-x1))+       d = (y1-y0)/(x1-x0)+       x = x0+   in  2*(dy0-d) / sqr (x0-x1) * (3*x-2*x1-x0) ++       2*(dy1-d) / sqr (x1-x0) * (3*x-2*x0-x1)
+ private/Numeric/Interpolation/Private/Piece.hs view
@@ -0,0 +1,31 @@+module Numeric.Interpolation.Private.Piece where++sqr :: (Num a) => a -> a+sqr x = x*x+++type T x y ny = (x, ny) -> (x, ny) -> x -> y++linear :: (Fractional a) => T a a a+linear (x0,y0) (x1,y1) x =+   (y0*(x1-x) + y1*(x-x0)) / (x1-x0)++hermite1 :: (Fractional a) => T a a (a, a)+hermite1 (x0,(y0,dy0)) (x1,(y1,dy1)) x =+   let d = (y1-y0)/(x1-x0)+   in  linear (x0,y0) (x1,y1) x ++       (dy0-d) * sqr ((x-x1)/(x0-x1)) * (x-x0) ++       (dy1-d) * sqr ((x-x0)/(x1-x0)) * (x-x1)++hermite1' :: (Fractional a) => T a a (a, a)+hermite1' (x0,(y0,dy0)) (x1,(y1,dy1)) x =+   let d = (y1-y0)/(x1-x0)+   in  d ++       (dy0-d) / sqr (x0-x1) * (2*(x-x1) * (x-x0) + sqr (x-x1)) ++       (dy1-d) / sqr (x1-x0) * (2*(x-x0) * (x-x1) + sqr (x-x0))++hermite1'' :: (Fractional a) => T a a (a, a)+hermite1'' (x0,(y0,dy0)) (x1,(y1,dy1)) x =+   let d = (y1-y0)/(x1-x0)+   in  2*(dy0-d) / sqr (x0-x1) * (3*x-2*x1-x0) ++       2*(dy1-d) / sqr (x1-x0) * (3*x-2*x0-x1)
+ src/Numeric/Interpolation/Basis.hs view
@@ -0,0 +1,62 @@+module Numeric.Interpolation.Basis (+   Compact.linear,+   Compact.hermite1,+   Compact.cubicLinear,+   Compact.cubicParabola,+   coefficientsToLinear,+   coefficientsToHermite1,+   coefficientsToCubicLinear,+   coefficientsToCubicParabola,+   ) where++import qualified Numeric.Interpolation.Basis.Compact as Compact+import qualified Numeric.Interpolation.NodeList as Nodes+import Numeric.Interpolation.Private.Basis+          (parabolaDerivativeCenterNode, hermite1Split)++import qualified Data.List as List+++{- |+@coefficientsToLinear nodes coefficients@+creates an interpolation function for @nodes@,+where the @coefficients@ correspond to the basis functions+constructed with @Basis.linear nodes@.+-}+coefficientsToLinear :: [a] -> [b] -> Nodes.T a b+coefficientsToLinear xs = Nodes.fromList . zip xs++{- |+Cf. 'coefficientsToLinear'+-}+coefficientsToHermite1 :: [a] -> [b] -> Nodes.T a (b, b)+coefficientsToHermite1 xs =+   Nodes.fromList . zip xs . hermite1Split xs++++mapAdjacent3 :: (a -> a -> a -> b) -> [a] -> [b]+mapAdjacent3 f xs0 =+   let xs1 = drop 1 xs0+       xs2 = drop 1 xs1+   in  List.zipWith3 f xs0 xs1 xs2++{- |+Cf. 'coefficientsToLinear'+-}+coefficientsToCubicLinear :: (Fractional a) => [a] -> [a] -> Nodes.T a (a, a)+coefficientsToCubicLinear xs =+   Nodes.fromList .+   mapAdjacent3 (\(xl,yl) (xn,yn) (xr,yr) -> (xn, (yn, (yr-yl)/(xr-xl)))) .+   zip xs++{- |+Cf. 'coefficientsToLinear'+-}+coefficientsToCubicParabola :: (Fractional a) => [a] -> [a] -> Nodes.T a (a, a)+coefficientsToCubicParabola xs =+   Nodes.fromList .+   mapAdjacent3+      (\pl pn@(xn,yn) pr ->+         (xn, (yn, parabolaDerivativeCenterNode pl pn pr))) .+   zip xs
+ src/Numeric/Interpolation/Basis/Compact.hs view
@@ -0,0 +1,109 @@+module Numeric.Interpolation.Basis.Compact (+   linear, hermite1, cubicLinear, cubicParabola,+   ) where++import qualified Numeric.Interpolation.NodeList as Nodes+import Numeric.Interpolation.Private.Basis (+   parabolaBasisDerivativeRight,+   parabolaBasisDerivativeCenter,+   parabolaBasisDerivativeLeft,+   )++import Control.Monad (liftM, liftM2)++import qualified Data.List as List+import Data.Maybe (catMaybes)+++mapAdjacentMaybe3 :: (Maybe a -> a -> Maybe a -> b) -> [a] -> [b]+mapAdjacentMaybe3 f xs =+   let jxs = map Just xs+   in  zipWith3 f (Nothing : jxs) xs (drop 1 jxs ++ [Nothing])++generic :: ny -> ny -> [x] -> [Nodes.T x ny]+generic nz ny =+   mapAdjacentMaybe3+      (\l n r ->+          Nodes.Node (n,ny)+             (maybe Nodes.Interval (flip Nodes.singleton nz) l)+             (maybe Nodes.Interval (flip Nodes.singleton nz) r))+++linear :: (Num b) => [a] -> [Nodes.T a b]+linear = generic 0 1++hermite1 :: (Num b) => [a] -> [Nodes.T a (b, b)]+hermite1 xs =+   generic (0,0) (1,0) xs+   +++   generic (0,0) (0,1) xs+++++mapAdjacentMaybe5 ::+   (Maybe a -> Maybe a -> a -> Maybe a -> Maybe a -> b) ->+   [a] -> [b]+mapAdjacentMaybe5 f xs =+   let jxs = map Just xs+       lxs1 = Nothing : jxs+       lxs2 = Nothing : lxs1+       rxs1 = drop 1 $ jxs ++ repeat Nothing+       rxs2 = drop 1 $ rxs1+   in  List.zipWith5 f lxs2 lxs1 xs rxs1 rxs2++cubicAutoGeneric ::+   (Num b) =>+   (a -> a -> a -> b) ->+   (a -> a -> a -> b) ->+   (a -> a -> a -> b) ->+   [a] -> [Nodes.T a (b, b)]+cubicAutoGeneric dl dn dr =+   mapAdjacentMaybe5+      (\ml2 ml1 n mr1 mr2 ->+         let node x y y' = (x, (y,y'))+         in  Nodes.fromList $ catMaybes $+             liftM (\l2 -> node l2 0 0) ml2 :+             liftM2 (\l2 l1 -> node l1 0 (dl l2 l1 n)) ml2 ml1 :+             liftM2 (\l1 r1 -> node n  1 (dn l1 n r1)) ml1 mr1 :+             liftM2 (\r1 r2 -> node r1 0 (dr n r1 r2)) mr1 mr2 :+             liftM (\r2 -> node r2 0 0) mr2 :+             [])+++{- |+Cubic interpolation+where the derivative at a node is set to the slope of the two adjacent nodes.+-}+cubicLinear :: (Fractional a) => [a] -> [Nodes.T a (a, a)]+cubicLinear =+   cubicAutoGeneric+      (\ll _l n -> recip $ n-ll)+      (\_l _n _r -> 0)+      (\n _r rr -> recip $ n-rr)+++{- |+Cubic interpolation+where the derivative at a node is set to the slope of the parabola+through the current and the two adjacent nodes.+-}+cubicParabola :: (Fractional a) => [a] -> [Nodes.T a (a, a)]+cubicParabola =+   cubicAutoGeneric+      parabolaBasisDerivativeRight+      parabolaBasisDerivativeCenter+      parabolaBasisDerivativeLeft+++{- |+Experimental interpolation+which is mean of 'cubicLinear' and 'cubicParabola'.+The result looks reasonable, too.+-}+_cubicMean :: (Fractional a) => [a] -> [Nodes.T a (a, a)]+_cubicMean =+   cubicAutoGeneric+      (\ll l n -> (parabolaBasisDerivativeRight ll l n + recip (n-ll))/2)+      (\l n r -> parabolaBasisDerivativeCenter l n r / 2)+      (\n r rr -> (parabolaBasisDerivativeLeft n r rr + recip (n-rr))/2)
+ src/Numeric/Interpolation/Basis/Full.hs view
@@ -0,0 +1,20 @@+module Numeric.Interpolation.Basis.Full (linear, hermite1) where++import qualified Numeric.Interpolation.NodeList as Nodes++import qualified Data.List.Match as Match+++generic :: ny -> ny -> [x] -> [Nodes.T x ny]+generic nz ny xs =+   map (Nodes.fromList . zip xs) $+   Match.take xs $ iterate (nz:) $ ny : repeat nz++linear :: (Num b) => [a] -> [Nodes.T a b]+linear = generic 0 1++hermite1 :: (Num b) => [a] -> [Nodes.T a (b, b)]+hermite1 xs =+   generic (0,0) (1,0) xs+   +++   generic (0,0) (0,1) xs
+ src/Numeric/Interpolation/NodeList.hs view
@@ -0,0 +1,48 @@+module Numeric.Interpolation.NodeList (+   T(Interval, Node),+   fromList,+   toList,+   singleton,+   lookup,+   ) where++import Data.Tuple.HT (mapFst)++import Prelude hiding (lookup)+++data T x y = Interval | Node (x, y) (T x y) (T x y)+   deriving (Eq, Ord, Show)++{- |+list must be sorted with respect to first element+-}+fromList :: [(x,y)] -> T x y+fromList =+   let merge n0 xys0 =+          case xys0 of+             (xy0,n1):(xy1,n2):xys ->+                (Node xy0 n0 n1,+                 uncurry (:) $ mapFst ((,) xy1) $ merge n2 xys)+             (xy0,n1):[] -> (Node xy0 n0 n1, [])+             [] -> (n0, [])+       rep (n,xyns) = if null xyns then n else rep $ merge n xyns+   in  rep . merge Interval . map (flip (,) Interval)++singleton :: x -> y -> T x y+singleton x y = Node (x,y) Interval Interval++toList :: T x y -> [(x,y)]+toList =+   let go Interval = []+       go (Node p l r) = go l ++ p : go r+   in  go++lookup :: Ord x => T x y -> x -> (Maybe (x,y), Maybe (x,y))+lookup nodes0 x0 =+   let go lb rb Interval = (lb, rb)+       go lb rb (Node n@(x,_y) ln rn) =+          if x0>=x+            then go (Just n) rb rn+            else go lb (Just n) ln+   in  go Nothing Nothing nodes0
+ src/Numeric/Interpolation/Piece.hs view
@@ -0,0 +1,22 @@+module Numeric.Interpolation.Piece (+   Piece.T,+   Piece.linear,+   hermite1,+   ) where++import qualified Numeric.Interpolation.Private.Piece as Piece+++{- |+Hermite interpolation with one derivative per node.+That is, the interpolating polynomial is cubic.+-}+hermite1 :: (Fractional a) => Piece.T a a (a, a)+hermite1 (x0,(y0,dy0)) (x1,(y1,dy1)) x =+   let d = (y1-y0) / dx10+       dx0 = x-x0+       dx1 = x1-x+       dx10 = x1-x0+   in  (y0*dx1 + y1*dx0 ++        ((dy0-d) * dx1 - (dy1-d) * dx0) * dx0 * dx1 / dx10)+          / dx10
+ src/Numeric/Interpolation/Piecewise.hs view
@@ -0,0 +1,30 @@+module Numeric.Interpolation.Piecewise (+   interpolate,+   interpolateConstantExt,+   ) where++import qualified Numeric.Interpolation.NodeList as Nodes+import qualified Numeric.Interpolation.Type as Type+++{- |+It is a checked error to interpolate outside of the range of nodes.+-}+interpolate :: (Ord x) => Type.T x y ny -> Nodes.T x ny -> x -> y+interpolate typ ns x =+   case Nodes.lookup ns x of+      (Just p0, Just p1) -> Type.interpolatePiece typ p0 p1 x+      _ -> error "interpolate: argument outside range"++{- |+Outside the range of nodes the interpolation function+takes the value of the respective border.+-}+interpolateConstantExt ::+   (Ord x) => Type.T x y ny -> Nodes.T x ny -> x -> y+interpolateConstantExt typ ns x =+   case Nodes.lookup ns x of+      (Just p0, Just p1) -> Type.interpolatePiece typ p0 p1 x+      (Just p, Nothing) -> Type.valueFromNode typ $ snd p+      (Nothing, Just p) -> Type.valueFromNode typ $ snd p+      (Nothing, Nothing) -> error "interpolateConstantExt: empty node list"
+ src/Numeric/Interpolation/Type.hs view
@@ -0,0 +1,81 @@+module Numeric.Interpolation.Type (+   T(..),+   linear,+   cubic,+   cubicLinear,+   cubicParabola,+   ) where++import qualified Numeric.Interpolation.NodeList as Nodes+import qualified Numeric.Interpolation.Piece as Piece+import qualified Numeric.Interpolation.Basis as Basis+import Numeric.Interpolation.Private.Basis (hermite1Split)+++data T x y ny =+   Cons {+      ssvFromNodes :: [x] -> [y] -> String,+      interpolatePiece :: Piece.T x y ny,+      basisFunctions :: [x] -> [Nodes.T x ny],+      coefficientsToInterpolator :: [x] -> [y] -> Nodes.T x ny,+      valueFromNode :: ny -> y+   }++linear :: T Double Double Double+linear =+   Cons {+      ssvFromNodes =+         \xs ys -> unlines $ zipWith (\x y -> show x ++ " " ++ show y) xs ys,+      interpolatePiece = Piece.linear,+      basisFunctions = Basis.linear,+      coefficientsToInterpolator = Basis.coefficientsToLinear,+      valueFromNode = id+   }++cubic :: T Double Double (Double, Double)+cubic =+   Cons {+      ssvFromNodes =+         \xs ys ->+            unlines .+            zipWith (\x (y,dy) -> show x ++ " " ++ show y ++ " " ++ show dy) xs $+            hermite1Split xs ys,+      interpolatePiece = Piece.hermite1,+      basisFunctions = Basis.hermite1,+      coefficientsToInterpolator = Basis.coefficientsToHermite1,+      valueFromNode = fst+   }++cubicLinear :: T Double Double (Double, Double)+cubicLinear =+   Cons {+      ssvFromNodes =+         \xs ys -> unlines $ zipWith (\x y -> show x ++ " " ++ show y) xs ys,+      interpolatePiece = Piece.hermite1,+      basisFunctions = Basis.cubicLinear,+      coefficientsToInterpolator = Basis.coefficientsToCubicLinear,+      valueFromNode = fst+   }++cubicParabola :: T Double Double (Double, Double)+cubicParabola =+   Cons {+      ssvFromNodes =+         \xs ys -> unlines $ zipWith (\x y -> show x ++ " " ++ show y) xs ys,+      interpolatePiece = Piece.hermite1,+      basisFunctions = Basis.cubicParabola,+      coefficientsToInterpolator = Basis.coefficientsToCubicParabola,+      valueFromNode = fst+   }+++_cubicMean :: T Double Double (Double, Double)+_cubicMean =+   Cons {+      ssvFromNodes =+         \xs ys -> unlines $ zipWith (\x y -> show x ++ " " ++ show y) xs ys,+      interpolatePiece = Piece.hermite1,+      basisFunctions = Basis.cubicParabola, -- Basis.cubicMean,+      coefficientsToInterpolator = Basis.coefficientsToCubicParabola, -- not correct+      valueFromNode = fst+   }
+ test/Test.hs view
@@ -0,0 +1,12 @@+module Main where++import qualified Test.Piece as Piece+++run :: String -> [(String, IO ())] -> IO ()+run prefix =+   mapM_ (\(msg,act) -> putStr (prefix ++ '.' : msg ++ ": ") >> act)++main :: IO ()+main =+   run "Piece" Piece.tests
+ test/Test/Piece.hs view
@@ -0,0 +1,60 @@+module Test.Piece where++import qualified Numeric.Interpolation.Piece as Piece+import qualified Numeric.Interpolation.Private.Piece as PiecePriv++import Test.QuickCheck (Property, quickCheck, (==>), )+++type Point = (Rational, Rational)++linearCommutative ::+   Point -> Point -> Rational -> Property+linearCommutative p1@(x1,_) p2@(x2,_) x =+   x1/=x2+   ==>+   Piece.linear p1 p2 x+   ==+   Piece.linear p2 p1 x+++type PointSlope = (Rational, (Rational, Rational))++hermite1Commutative ::+   PointSlope -> PointSlope -> Rational -> Property+hermite1Commutative p1@(x1,_) p2@(x2,_) x =+   x1/=x2+   ==>+   Piece.hermite1 p1 p2 x+   ==+   Piece.hermite1 p2 p1 x+++linearHermite1 ::+   Point -> Point -> Rational -> Property+linearHermite1 p1@(x1,y1) p2@(x2,y2) x =+   x1/=x2+   ==>+   Piece.linear p1 p2 x+   ==+   let slope = (y2-y1)/(x2-x1)+   in  Piece.hermite1 (x1, (y1,slope)) (x2, (y2, slope)) x+++hermite1Alternative ::+   PointSlope -> PointSlope -> Rational -> Property+hermite1Alternative p1@(x1,_) p2@(x2,_) x =+   x1/=x2+   ==>+   Piece.hermite1 p1 p2 x+   ==+   PiecePriv.hermite1 p1 p2 x+++tests :: [(String, IO ())]+tests =+   ("linearCommutative", quickCheck linearCommutative) :+   ("hermite1Commutative", quickCheck hermite1Commutative) :+   ("linearHermite1", quickCheck linearHermite1) :+   ("hermite1Alternative", quickCheck hermite1Alternative) :+   []