diff --git a/Control/Monad/Numeric.hs b/Control/Monad/Numeric.hs new file mode 100644 --- /dev/null +++ b/Control/Monad/Numeric.hs @@ -0,0 +1,42 @@ +-- | +-- Module : Control.Monad.Numeric +-- Copyright : (c) 2011 Aleksey Khudyakov +-- License : BSD3 +-- +-- Maintainer : Aleksey Khudyakov +-- 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 #-} diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Numeric/Classes/Indexing.hs b/Numeric/Classes/Indexing.hs new file mode 100644 --- /dev/null +++ b/Numeric/Classes/Indexing.hs @@ -0,0 +1,56 @@ +{-# LANGUAGE TypeFamilies #-} +-- | +-- Module : Numeric.Classes.Indexing +-- Copyright : (c) 2011 Aleksey Khudyakov +-- License : BSD3 +-- +-- Maintainer : Aleksey Khudyakov +-- 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.!) diff --git a/Numeric/Tools/Differentiation.hs b/Numeric/Tools/Differentiation.hs new file mode 100644 --- /dev/null +++ b/Numeric/Tools/Differentiation.hs @@ -0,0 +1,146 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE ForeignFunctionInterface #-} +-- | +-- Module : Numeric.Tools.Differentiation +-- Copyright : (c) 2011 Aleksey Khudyakov +-- License : BSD3 +-- +-- Maintainer : Aleksey Khudyakov +-- 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. diff --git a/Numeric/Tools/Equation.hs b/Numeric/Tools/Equation.hs new file mode 100644 --- /dev/null +++ b/Numeric/Tools/Equation.hs @@ -0,0 +1,44 @@ +-- | +-- Module : Numeric.Tools.Equation +-- Copyright : (c) 2011 Aleksey Khudyakov +-- License : BSD3 +-- +-- Maintainer : Aleksey Khudyakov +-- 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 + diff --git a/Numeric/Tools/Integration.hs b/Numeric/Tools/Integration.hs new file mode 100644 --- /dev/null +++ b/Numeric/Tools/Integration.hs @@ -0,0 +1,184 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE DeriveDataTypeable #-} +-- | +-- Module : Numeric.Tools.Integration +-- Copyright : (c) 2011 Aleksey Khudyakov +-- License : BSD3 +-- +-- Maintainer : Aleksey Khudyakov +-- 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 diff --git a/Numeric/Tools/Interpolation.hs b/Numeric/Tools/Interpolation.hs new file mode 100644 --- /dev/null +++ b/Numeric/Tools/Interpolation.hs @@ -0,0 +1,197 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE TypeFamilies #-} +-- | +-- Module : Numeric.Tools.Interpolation +-- Copyright : (c) 2011 Aleksey Khudyakov +-- License : BSD3 +-- +-- Maintainer : Aleksey Khudyakov +-- 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 #-} diff --git a/Numeric/Tools/Mesh.hs b/Numeric/Tools/Mesh.hs new file mode 100644 --- /dev/null +++ b/Numeric/Tools/Mesh.hs @@ -0,0 +1,81 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE TypeFamilies #-} +-- | +-- Module : Numeric.Tools.Mesh +-- Copyright : (c) 2011 Aleksey Khudyakov +-- License : BSD3 +-- +-- Maintainer : Aleksey Khudyakov +-- 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 diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/cbits/ieee.c b/cbits/ieee.c new file mode 100644 --- /dev/null +++ b/cbits/ieee.c @@ -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; +} diff --git a/numeric-tools.cabal b/numeric-tools.cabal new file mode 100644 --- /dev/null +++ b/numeric-tools.cabal @@ -0,0 +1,34 @@ +Name: numeric-tools +Version: 0.1.0.0 +Cabal-Version: >= 1.6 +License: BSD3 +License-File: LICENSE +Author: Aleksey Khudyakov +Maintainer: Aleksey Khudyakov +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