numeric-tools (empty) → 0.1.0.0
raw patch · 11 files changed
+823/−0 lines, 11 filesdep +basedep +ieee754dep +vectorsetup-changed
Dependencies added: base, ieee754, vector
Files
- Control/Monad/Numeric.hs +42/−0
- LICENSE +30/−0
- Numeric/Classes/Indexing.hs +56/−0
- Numeric/Tools/Differentiation.hs +146/−0
- Numeric/Tools/Equation.hs +44/−0
- Numeric/Tools/Integration.hs +184/−0
- Numeric/Tools/Interpolation.hs +197/−0
- Numeric/Tools/Mesh.hs +81/−0
- Setup.hs +2/−0
- cbits/ieee.c +7/−0
- numeric-tools.cabal +34/−0
+ Control/Monad/Numeric.hs view
@@ -0,0 +1,42 @@+-- |+-- Module : Control.Monad.Numeric+-- Copyright : (c) 2011 Aleksey Khudyakov+-- License : BSD3+--+-- Maintainer : Aleksey Khudyakov <alexey.skladnoy@gmail.com>+-- Stability : experimental+-- Portability : portable+--+-- Function useful for writing numeric code which works with mutable+-- data.+module Control.Monad.Numeric (+ forGen+ , for+ ) where++-- | For function which act much like for loop in the C+forGen :: Monad m + => a -- ^ Staring index value+ -> (a -> Bool) -- ^ Condition+ -> (a -> a) -- ^ Function to modify index+ -> (a -> m ()) -- ^ Action to perform+ -> m ()+forGen n test next a = worker n+ where+ worker i | test i = a i >> worker (next i)+ | otherwise = return ()+{-# INLINE forGen #-}++-- | Specialized for loop. Akin to:+--+-- > for( i = 0; i < 10; i++) { ...+for :: Monad m + => Int -- ^ Starting index+ -> Int -- ^ Maximal index value not reached+ -> (Int -> m ()) -- ^ Action to perfor,+ -> m ()+for i maxI a = worker i+ where+ worker j | j < maxI = a j >> worker (j+1)+ | otherwise = return ()+{-# INLINE for #-}
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) Alexey Khudyakov++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 author nor the names of his contributors+ may be used to endorse or promote products derived from this software+ without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE 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 AUTHORS 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.
+ Numeric/Classes/Indexing.hs view
@@ -0,0 +1,56 @@+{-# LANGUAGE TypeFamilies #-}+-- |+-- Module : Numeric.Classes.Indexing+-- Copyright : (c) 2011 Aleksey Khudyakov+-- License : BSD3+--+-- Maintainer : Aleksey Khudyakov <alexey.skladnoy@gmail.com>+-- Stability : experimental+-- Portability : portable+--+module Numeric.Classes.Indexing (+ Indexable(..)+ , validIndex+ ) where++import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Storable as S++++-- | Type class for array-like data type which support @O(1)@ access+-- by integer index starting from zero.+class Indexable a where+ type IndexVal a :: *+ -- | Size of table.+ size :: a -> Int+ -- | /O(1)/ Index table without range cheking.+ unsafeIndex :: a -> Int -> IndexVal a+ -- | /O(1)/ Safe indexing. Calls error if index is out of range.+ (!) :: a -> Int -> IndexVal a+ x ! i | i < 0 || i > size x = error "Numeric.Classes.Indexing.!: index is out of range"+ | otherwise = unsafeIndex x i++-- | Check that index is valid+validIndex :: Indexable a => a -> Int -> Bool +validIndex tbl i = i >= 0 && i < size tbl+{-# INLINE validIndex #-}++instance Indexable (V.Vector a) where+ type IndexVal (V.Vector a) = a+ size = V.length+ unsafeIndex = V.unsafeIndex+ (!) = (V.!)++instance U.Unbox a => Indexable (U.Vector a) where+ type IndexVal (U.Vector a) = a+ size = U.length+ unsafeIndex = U.unsafeIndex+ (!) = (U.!)++instance S.Storable a => Indexable (S.Vector a) where+ type IndexVal (S.Vector a) = a+ size = S.length+ unsafeIndex = S.unsafeIndex+ (!) = (S.!)
+ Numeric/Tools/Differentiation.hs view
@@ -0,0 +1,146 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE ForeignFunctionInterface #-}+-- |+-- Module : Numeric.Tools.Differentiation+-- Copyright : (c) 2011 Aleksey Khudyakov+-- License : BSD3+--+-- Maintainer : Aleksey Khudyakov <alexey.skladnoy@gmail.com>+-- Stability : experimental+-- Portability : portable+--+-- Numerical differentiation. 'diffRichardson' is preferred way to+-- calculate derivative.+--+module Numeric.Tools.Differentiation (+ -- * Differentiation+ DiffRes(..)+ , diffRichardson+ -- * Fast but imprecise+ , diffSimple+ , diffSimmetric+ -- * Utils+ , representableDelta + -- * References+ -- $references+ ) where++import Control.Monad.ST (runST)+import Data.Data (Data,Typeable)+import qualified Data.Vector.Unboxed.Mutable as M+import Foreign+import Foreign.C++import Numeric.IEEE (infinity, nan)++++-- | Differentiation result+data DiffRes = DiffRes { diffRes :: Double -- ^ Derivative value+ , diffPrecision :: Double -- ^ Rough error estimate+ }+ deriving (Show,Eq,Data,Typeable)++-- | Calculate derivative using Richaradson's deferred approach to+-- limit. This is a preferred method for numeric differentiation+-- since it's most precise. Function could be evaluated up to 20+-- times.+--+-- Initial step size should be chosen fairly big. Too small one will+-- result reduced precision, too big one in nonsensical answer.+diffRichardson :: (Double -> Double) -- ^ Function+ -> Double -- ^ Delta+ -> Double -- ^ Point at which evaluate differential+ -> DiffRes+diffRichardson f h x0 = runST $ do+ let nMax = 10 -- Maximum number of iterations+ let con = 1.4 -- Decrement for step size+ con2 = con*con -- Square of decrement+ let safe = 2+ -- Start calculations+ arr <- M.new nMax+ let worker i hh err ans = do+ -- Calculate extrapolations+ let richard j fac x err' ans' = do+ xOld <- replace arr (j-1) x+ case () of+ _| j > i -> return (ans',err')+ | otherwise -> + let x' = (x*fac - xOld) / (fac - 1) -- New extrapolation+ errt = max (abs $ x' - x) (abs $ x' - xOld) -- New error estimate+ (ans'',err'') = if errt < err' then (x' , errt)+ else (ans' , err')+ in richard (j+1) (fac*con2) x' err'' ans''+ -- Main loop+ let hh' = hh / con -- New step size+ d = (f (x0 + hh') - f (x0 - hh')) / (2 * hh') -- New approximation+ x' <- M.read arr (i-1)+ (ans',err') <- richard 1 con2 d err ans+ x'' <- M.read arr i+ case () of+ _| abs (x' - x'') >= safe * err' -> return $ DiffRes ans' err'+ | i >= nMax - 1 -> return $ DiffRes ans' err'+ | otherwise -> worker (i+1) hh' err' ans'+ -- Calculate+ M.write arr 0 $ (f (x0 + h) - f (x0 - h)) / (2*h)+ worker 1 h infinity nan++++-- | Simplest form of differentiation. Should be used only when+-- function evaluation is prohibitively expensive and already+-- computed value at point @x@ should be reused.+--+-- > f'(x) = f(x+h) - f(x) / h+diffSimple :: (Double -> Double) -- ^ Function to differentiate+ -> Double -- ^ Delta+ -> (Double,Double) -- ^ Coordinate and function value at this point+ -> Double+diffSimple f h (x,fx) = (f (x + h') - fx) / h' where h' = representableDelta x h+{-# INLINE diffSimple #-} +++-- | Simple differentiation. It uses simmetric rule and provide+-- reasonable accuracy. It's suitable when function evaluation is+-- expensive and precision could be traded for speed.+--+-- > f'(x) = f(x-h) + f(x+h) / 2h+diffSimmetric :: (Double -> Double) -- ^ Function to differentiate+ -> Double -- ^ Delta+ -> Double -- ^ Point at which evaluate differential+ -> Double+diffSimmetric f h x = (f(x + h') - f(x - h')) / (2 * h')+ where+ h' = representableDelta x h+++ +----------------------------------------------------------------+-- Helpers+----------------------------------------------------------------++-- replace :: (PrimMonad m, M.MVector v a) => v (PrimState m) a -> Int -> a -> m a+replace arr i x = do+ x' <- M.read arr i+ M.write arr i x+ return x'+{-# INLINE replace #-}+ ++-- | For number @x@ and small @h@ return such @h'@ that @x+h'@ and @x@+-- differ by representable number+representableDelta :: Double -- ^ x+ -> Double -- ^ small delta+ -> Double +representableDelta x h = realToFrac $ unsafePerformIO $ representableDeltaFFI (realToFrac x) (realToFrac h)+{-# INLINE representableDelta #-}++foreign import ccall "numeric_tools_representable_delta" + representableDeltaFFI :: CDouble -> CDouble -> IO CDouble+++-- $references+--+-- * Ridders, C.J.F. 1982, Accurate computation of F`(x) and+-- F`(x)F``(x), Advances in Engineering Software, vol. 4, no. 2,+-- pp. 75-76.
+ Numeric/Tools/Equation.hs view
@@ -0,0 +1,44 @@+-- |+-- Module : Numeric.Tools.Equation+-- Copyright : (c) 2011 Aleksey Khudyakov+-- License : BSD3+--+-- Maintainer : Aleksey Khudyakov <alexey.skladnoy@gmail.com>+-- Stability : experimental+-- Portability : portable+--+-- Numerical solution of ordinary equations.+module Numeric.Tools.Equation ( + solveBisection+ ) where++import Numeric.IEEE (epsilon)++++-- | Solve equation @f(x) = 0@ using bisection method. Function is+-- must be continous. If function has different signs at the ends of+-- initial interval answer is always returned. 'Nothing' is returned+-- if function fails to find an answer.+solveBisection :: Double -- ^ Required absolute precision+ -> (Double,Double) -- ^ Range+ -> (Double -> Double) -- ^ Equation+ -> Maybe Double+solveBisection eps (a,b) f+ | a >= b = Nothing+ | fa * fb > 0 = Nothing+ | otherwise = Just $ bisectionWorker (abs eps) f a b fa fb+ where+ fa = f a+ fb = f b++bisectionWorker :: Double -> (Double -> Double) -> Double -> Double -> Double -> Double -> Double+bisectionWorker eps f a b fa fb+ | (b - a) <= eps = c+ | (b - a) / b <= epsilon = c+ | fa * fc < 0 = bisectionWorker eps f a c fa fc+ | otherwise = bisectionWorker eps f c b fc fb+ where+ c = 0.5 * (a + b)+ fc = f c+
+ Numeric/Tools/Integration.hs view
@@ -0,0 +1,184 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DeriveDataTypeable #-}+-- |+-- Module : Numeric.Tools.Integration+-- Copyright : (c) 2011 Aleksey Khudyakov+-- License : BSD3+--+-- Maintainer : Aleksey Khudyakov <alexey.skladnoy@gmail.com>+-- Stability : experimental+-- Portability : portable+--+-- Funtions for numerical integration. 'quadRomberg' or 'quadSimpson'+-- are reasonable choices in most cases. For non-smooth function they+-- converge poorly and 'quadTrapezoid' should be used then.+--+-- For example this code intergrates exponent from 0 to 1:+--+-- >>> let res = quadRomberg defQuad (0, 1) exp+--+-- >>> quadRes res -- Integration result+-- Just 1.718281828459045+--+-- >>> quadPrecEst res -- Estimate of precision+-- 2.5844957590474064e-16+--+-- >>> quadNIter res -- Number of iterations performed+-- 6+module Numeric.Tools.Integration (+ -- * Integration parameters and results+ QuadParam(..)+ , defQuad+ , QuadRes(..)+ -- * Integration functions+ , quadTrapezoid+ , quadSimpson+ , quadRomberg+ ) where++import Control.Monad.ST++import Data.Data (Data,Typeable)+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as M++++----------------------------------------------------------------+-- Data types+----------------------------------------------------------------++-- | Integration parameters for numerical routines. Note that each+-- additional iteration doubles number of function evaluation required+-- to compute integral.+--+-- Number of iterations is capped at 30.+data QuadParam = QuadParam {+ quadPrecision :: Double -- ^ Relative precision of answer+ , quadMaxIter :: Int -- ^ Maximum number of iterations+ }+ deriving (Show,Eq,Data,Typeable)++-- Number of iterations limited to 30+maxIter :: QuadParam -> Int+maxIter = min 30 . quadMaxIter++-- | Default parameters for integration functions+--+-- * Maximum number of iterations = 20+--+-- * Precision is 10⁻⁹+defQuad :: QuadParam+defQuad = QuadParam { quadPrecision = 1e-9+ , quadMaxIter = 20+ }++-- | Result of numeric integration.+data QuadRes = QuadRes { quadRes :: Maybe Double -- ^ Integraion result+ , quadPrecEst :: Double -- ^ Rough estimate of attained precision+ , quadNIter :: Int -- ^ Number of iterations+ }+ deriving (Show,Eq,Data,Typeable)++++----------------------------------------------------------------+-- Different integration methods+----------------------------------------------------------------++-- | Integration of using trapezoids. This is robust algorithm and+-- place and useful for not very smooth. But it is very slow. It+-- hundreds times slower then 'quadRomberg' if function is+-- sufficiently smooth.+quadTrapezoid :: QuadParam -- ^ Parameters+ -> (Double, Double) -- ^ Integration limits+ -> (Double -> Double) -- ^ Function to integrate+ -> QuadRes+quadTrapezoid param (a,b) f = worker 1 1 (trapGuess a b f)+ where+ eps = quadPrecision param -- Requred precision+ maxN = maxIter param -- Maximum allowed number of iterations+ worker n nPoints q+ | n > 5 && d < eps = ret (Just q')+ | n >= maxN = ret Nothing+ | otherwise = worker (n+1) (nPoints*2) q'+ where+ q' = nextTrapezoid a b nPoints f q -- New approximation+ d = abs (q' - q) / abs q -- Precision estimate+ ret = \x -> QuadRes x d n++-- | Integration using Simpson rule. It should be more efficient than+-- 'quadTrapezoid' if function being integrated have finite fourth+-- derivative.+quadSimpson :: QuadParam -- ^ Parameters+ -> (Double, Double) -- ^ Integration limits+ -> (Double -> Double) -- ^ Function to integrate+ -> QuadRes+quadSimpson param (a,b) f = worker 1 1 0 (trapGuess a b f)+ where+ eps = quadPrecision param -- Requred precision+ maxN = maxIter param -- Maximum allowed number of points for evaluation+ worker n nPoints s st+ | n > 5 && d < eps = ret (Just s')+ | n >= maxN = ret Nothing+ | otherwise = worker (n+1) (nPoints*2) s' st'+ where+ st' = nextTrapezoid a b nPoints f st+ s' = (4*st' - st) / 3+ d = abs (s' - s) / abs s+ ret = \x -> QuadRes x d n++-- | Integration using Romberg rule. For sufficiently smooth functions+-- (e.g. analytic) it's a fastest of three.+quadRomberg :: QuadParam -- ^ Parameters+ -> (Double, Double) -- ^ Integration limits+ -> (Double -> Double) -- ^ Function to integrate+ -> QuadRes+quadRomberg param (a,b) f =+ runST $ do+ let eps = quadPrecision param+ maxN = maxIter param+ arr <- M.new maxN+ -- Calculate new approximation+ let nextAppr n = runNextAppr 0 4 where+ runNextAppr i fac s = do+ x <- M.read arr i+ M.write arr i s+ if i >= n+ then return s+ else runNextAppr (i+1) (fac*4) $ s + (s - x) / (fac - 1)+ -- Maine loop+ let worker n nPoints st s = do+ let st' = nextTrapezoid a b nPoints f st+ s' <- M.write arr 0 st >> nextAppr n st'+ let d = abs (s' - s) / abs s+ case () of+ _ | n > 5 && d < eps -> return $ QuadRes (Just s') d n+ | n >= maxN -> return $ QuadRes Nothing d n+ | otherwise -> worker (n+1) (nPoints*2) st' s'+ -- Calculate integral+ worker 1 1 st0 st0 where st0 = trapGuess a b f++++----------------------------------------------------------------+-- Helpers+----------------------------------------------------------------++-- Initial guess for trapezoid rule+trapGuess :: Double -> Double -> (Double -> Double) -> Double+trapGuess !a !b f = 0.5 * (b - a) * (f b + f a)+++-- Refinement of guess using trapeziod algorithms+nextTrapezoid :: Double -- Lower integration limit+ -> Double -- Upper integration limit+ -> Int -- Number of additional points+ -> (Double -> Double) -- Function to integrate+ -> Double -- Approximation+ -> Double+nextTrapezoid !a !b !n f !q = 0.5 * (q + sep * s)+ where+ sep = (b - a) / fromIntegral n -- Separation between points+ x0 = a + 0.5 * sep -- Starting point+ s = U.sum $ U.map f $ U.iterateN n (+sep) x0 -- Sum of all points
+ Numeric/Tools/Interpolation.hs view
@@ -0,0 +1,197 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE TypeFamilies #-}+-- |+-- Module : Numeric.Tools.Interpolation+-- Copyright : (c) 2011 Aleksey Khudyakov+-- License : BSD3+--+-- Maintainer : Aleksey Khudyakov <alexey.skladnoy@gmail.com>+-- Stability : experimental+-- Portability : portable+--+-- Function interpolation.+--+-- Sine interpolation using cubic splines:+--+-- >>> let tbl = cubicSpline $ tabulateFun (uniformMesh (0,10) 100) sin+-- >>> tbl `at` 1.786+-- 0.9769239849844867+module Numeric.Tools.Interpolation (+ -- * Type class+ Interpolation(..)+ , tabulate+ -- * Linear interpolation+ , LinearInterp+ , linearInterp+ -- * Cubic splines+ , CubicSpline+ , cubicSpline+ --+ , module Numeric.Tools.Mesh+ ) where++import Control.Monad.ST (runST)+import Data.Data (Data,Typeable)++import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as M++import Control.Monad.Numeric+import Numeric.Classes.Indexing+import Numeric.Tools.Mesh++++----------------------------------------------------------------++-- | Interpolation for arbitraty meshes+class Interpolation a where+ -- | Interpolate function at some point. Function should not+ -- fail outside of mesh however it may and most likely will give+ -- nonsensical results+ at :: (IndexVal m ~ Double, Mesh m) => a m -> Double -> Double+ -- | Tabulate function+ tabulateFun :: (IndexVal m ~ Double, Mesh m) => m -> (Double -> Double) -> a m+ -- | Use table of already evaluated function and mesh. Sizes of mesh+ -- and table must coincide but it's not checked. Do not use this+ -- function use 'tabulate' instead.+ unsafeTabulate :: (IndexVal m ~ Double, Mesh m, G.Vector v Double) => m -> v Double -> a m+ -- | Get mesh.+ interpolationMesh :: a m -> m+ -- | Get table of function values + interpolationTable :: a m -> U.Vector Double+ ++-- | Use table of already evaluated function and mesh. Sizes of mesh+-- and table must coincide. +tabulate :: (Interpolation a, IndexVal m ~ Double, Mesh m, G.Vector v Double) => m -> v Double -> a m+tabulate mesh tbl+ | size mesh /= G.length tbl = error "Numeric.Tools.Interpolation.tabulate: size of vector and mesh do not match"+ | otherwise = unsafeTabulate mesh tbl+{-# INLINE tabulate #-}++----------------------------------------------------------------+-- Linear interpolation+----------------------------------------------------------------++-- | Data for linear interpolation+data LinearInterp a = LinearInterp { linearInterpMesh :: a+ , linearInterpTable :: U.Vector Double+ }+ deriving (Show,Eq,Data,Typeable)++-- | Function used to fix types+linearInterp :: LinearInterp a -> LinearInterp a+linearInterp = id++instance Mesh a => Indexable (LinearInterp a) where+ type IndexVal (LinearInterp a) = (IndexVal a, Double)+ size (LinearInterp _ vec) = size vec+ unsafeIndex (LinearInterp mesh vec) i = ( unsafeIndex mesh i+ , unsafeIndex vec i+ )+ {-# INLINE size #-}+ {-# INLINE unsafeIndex #-}++instance Interpolation LinearInterp where+ at = linearInterpolation+ tabulateFun mesh f = LinearInterp mesh (U.generate (size mesh) (f . unsafeIndex mesh))+ unsafeTabulate mesh tbl = LinearInterp mesh (G.convert tbl)+ interpolationMesh = linearInterpMesh+ interpolationTable = linearInterpTable++linearInterpolation :: (Mesh a, IndexVal a ~ Double) => LinearInterp a -> Double -> Double+linearInterpolation tbl@(LinearInterp mesh _) x = a + (x - xa) / (xb - xa) * (b - a)+ where+ i = safeFindIndex mesh x+ (xa,a) = unsafeIndex tbl i+ (xb,b) = unsafeIndex tbl (i+1)++++----------------------------------------------------------------+-- Cubic splines+----------------------------------------------------------------++-- | Natural cubic splines+data CubicSpline a = CubicSpline { cubicSplineMesh :: a+ , cubicSplineTable :: U.Vector Double+ , cubicSplineY2 :: U.Vector Double+ }+ deriving (Eq,Show,Data,Typeable)++-- | Function used to fix types+cubicSpline :: CubicSpline a -> CubicSpline a +cubicSpline = id++instance Interpolation CubicSpline where+ at (CubicSpline mesh ys y2) x = y+ where+ i = safeFindIndex mesh x+ -- Table lookup+ xa = unsafeIndex mesh i+ xb = unsafeIndex mesh (i+1)+ ya = unsafeIndex ys i+ yb = unsafeIndex ys (i+1)+ da = unsafeIndex y2 i+ db = unsafeIndex y2 (i+1)+ -- + h = xb - xa+ a = (xb - x ) / h+ b = (x - xa) / h+ y = a * ya + b * yb + + ((a*a*a - a) * da + (b*b*b - b) * db) * (h * h) / 6+ ------+ tabulateFun mesh f = makeCubicSpline mesh (U.generate (size mesh) (f . unsafeIndex mesh))+ unsafeTabulate mesh tbl = makeCubicSpline mesh (G.convert tbl)+ interpolationMesh = cubicSplineMesh+ interpolationTable = cubicSplineTable+ ++-- These are natural cubic splines+makeCubicSpline :: (IndexVal a ~ Double, Mesh a) => a -> U.Vector Double -> CubicSpline a+makeCubicSpline xs ys = runST $ do+ let n = size ys+ y2 <- M.new n+ u <- M.new n+ M.write y2 0 0.0+ M.write u 0 0.0+ -- Forward pass+ for 1 (n-1) $ \i -> do+ yVal <- M.read y2 (i-1)+ uVal <- M.read u (i-1)+ let sig = delta xs i / delta xs (i+1)+ p = sig * yVal + 2+ u' = delta ys (i+1) / delta xs (i+1) - delta ys i / delta xs i+ M.write y2 i $ (sig - 1) / p+ M.write u i $ (6 * u' / (xs ! (i+1) - xs ! (i-1)) - sig * uVal) / p+ -- Backward pass+ M.write y2 (n-1) 0.0+ forGen (n-2) (>= 0) pred $ \i -> do+ uVal <- M.read u i+ yVal <- M.read y2 i+ yVal1 <- M.read y2 (i+1)+ M.write y2 i $ yVal * yVal1 + uVal+ -- Done+ y2' <- G.unsafeFreeze y2+ return (CubicSpline xs ys y2')+++----------------------------------------------------------------+-- Helpers++delta :: (Num (IndexVal a), Indexable a) => a -> Int -> IndexVal a+delta tbl i = (tbl ! i) - (tbl ! (i - 1))+{-# INLINE delta #-}++safeFindIndex :: Mesh a => a -> Double -> Int+safeFindIndex mesh x = + case meshFindIndex mesh x of+ i | i < 0 -> 0+ | i > n -> n+ | otherwise -> i+ where+ n = size mesh - 2+{-# INLINE safeFindIndex #-}
+ Numeric/Tools/Mesh.hs view
@@ -0,0 +1,81 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE TypeFamilies #-}+-- |+-- Module : Numeric.Tools.Mesh+-- Copyright : (c) 2011 Aleksey Khudyakov+-- License : BSD3+--+-- Maintainer : Aleksey Khudyakov <alexey.skladnoy@gmail.com>+-- Stability : experimental+-- Portability : portable+--+-- 1-dimensional meshes. Used by 'Numeric.Tools.Interpolation'.+--+module Numeric.Tools.Mesh (+ -- * Meshes+ Mesh(..)+ -- ** Uniform mesh+ , UniformMesh+ , uniformMesh+ , uniformMeshStep+ ) where++import Data.Data (Data,Typeable)+import Numeric.Classes.Indexing++++----------------------------------------------------------------+-- Type class+----------------------------------------------------------------++-- | Class for 1-dimensional meshes. Mesh is ordered set of+-- points. Each instance must guarantee that every next point is+-- greater that previous and there is at least 2 points in mesh.+class Indexable a => Mesh a where+ -- | Low bound of mesh+ meshLowerBound :: a -> Double+ -- | Upper bound of mesh+ meshUpperBound :: a -> Double++ -- | Find such index for value that+ --+ -- > mesh ! i <= x && mesh ! i+1 > x+ --+ -- Will return invalid index if value is out of range.+ meshFindIndex :: a -> Double -> Int+++++----------------------------------------------------------------+-- Uniform mesh+----------------------------------------------------------------++-- | Uniform mesh+data UniformMesh = UniformMesh { uniformMeshFrom :: Double+ , uniformMeshStep :: Double + -- ^ Distance between points+ , uniformMeshSize :: Int+ }+ deriving (Eq,Show,Data,Typeable)++-- | Create uniform mesh+uniformMesh :: (Double,Double) -- ^ Lower and upper bound+ -> Int -- ^ Number of points+ -> UniformMesh+uniformMesh (a,b) n+ | b <= a = error "Numeric.Tool.Interpolation.Mesh.uniformMesh: bad range"+ | n < 2 = error "Numeric.Tool.Interpolation.Mesh.uniformMesh: too few points"+ | otherwise = UniformMesh a ((b - a) / fromIntegral (n - 1)) n+++instance Indexable UniformMesh where+ type IndexVal UniformMesh = Double+ size = uniformMeshSize+ unsafeIndex (UniformMesh a da _) i = a + fromIntegral i * da++instance Mesh UniformMesh where+ meshLowerBound = uniformMeshFrom+ meshUpperBound (UniformMesh a da n) = a + da * fromIntegral (n - 1)+ meshFindIndex (UniformMesh a da _) x = truncate $ (x - a) / da
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ cbits/ieee.c view
@@ -0,0 +1,7 @@++double numeric_tools_representable_delta(double x, double h)+{+ /* temp is volatile to force loading from registers to memory. */+ volatile double temp = x + h;+ return temp - x;+}
+ numeric-tools.cabal view
@@ -0,0 +1,34 @@+Name: numeric-tools+Version: 0.1.0.0+Cabal-Version: >= 1.6+License: BSD3+License-File: LICENSE+Author: Aleksey Khudyakov <alexey.skladnoy@gmail.com>+Maintainer: Aleksey Khudyakov <alexey.skladnoy@gmail.com>+Homepage: https://bitbucket.org/Shimuuar/numeric-tools+bug-reports: https://bitbucket.org/Shimuuar/numeric-tools/issues+Category: Math, Numerical+Build-Type: Simple+Synopsis: Collection of numerical tools for integration, differentiation etc.+ +Description:+ Package provides function to perform numeric integration and+ differentiation, function interpolation.++source-repository head+ type: hg+ location: https://bitbucket.org/Shimuuar/numeric-tools++Library+ Build-Depends: base >=3 && <5,+ ieee754 >= 0.7.3,+ vector >= 0.7.0.1+ Exposed-modules: Control.Monad.Numeric+ Numeric.Classes.Indexing+ Numeric.Tools.Equation+ Numeric.Tools.Differentiation+ Numeric.Tools.Integration+ Numeric.Tools.Interpolation+ Numeric.Tools.Mesh+ c-sources: cbits/ieee.c+ ghc-options: -Wall -O2