diff --git a/hmatrix-gsl.cabal b/hmatrix-gsl.cabal --- a/hmatrix-gsl.cabal +++ b/hmatrix-gsl.cabal @@ -1,5 +1,5 @@ Name: hmatrix-gsl -Version: 0.16.0.3 +Version: 0.17.0.0 License: GPL License-file: LICENSE Author: Alberto Ruiz @@ -25,7 +25,7 @@ library - Build-Depends: base<5, hmatrix>=0.16, array, vector, + Build-Depends: base<5, hmatrix>=0.17, array, vector, process, random @@ -43,6 +43,8 @@ Numeric.GSL.ODE, Numeric.GSL, Numeric.GSL.LinearAlgebra, + Numeric.GSL.Interpolation, + Numeric.GSL.SimulatedAnnealing, Graphics.Plot other-modules: Numeric.GSL.Internal, Numeric.GSL.Vector, @@ -52,7 +54,12 @@ C-sources: src/Numeric/GSL/gsl-aux.c - cc-options: -O4 -msse2 -Wall + cc-options: -O4 -Wall + + if arch(x86_64) + cc-options: -msse2 + if arch(i386) + cc-options: -msse2 ghc-options: -Wall -fno-warn-missing-signatures -fno-warn-orphans diff --git a/src/Graphics/Plot.hs b/src/Graphics/Plot.hs --- a/src/Graphics/Plot.hs +++ b/src/Graphics/Plot.hs @@ -27,13 +27,13 @@ ) where -import Numeric.Container +import Numeric.LinearAlgebra.HMatrix import Data.List(intersperse) import System.Process (system) -- | From vectors x and y, it generates a pair of matrices to be used as x and y arguments for matrix functions. meshdom :: Vector Double -> Vector Double -> (Matrix Double , Matrix Double) -meshdom r1 r2 = (outer r1 (constant 1 (dim r2)), outer (constant 1 (dim r1)) r2) +meshdom r1 r2 = (outer r1 (konst 1 (size r2)), outer (konst 1 (size r1)) r2) {- | Draws a 3D surface representation of a real matrix. diff --git a/src/Numeric/GSL.hs b/src/Numeric/GSL.hs --- a/src/Numeric/GSL.hs +++ b/src/Numeric/GSL.hs @@ -22,6 +22,7 @@ , module Numeric.GSL.Root , module Numeric.GSL.ODE , module Numeric.GSL.Fitting +, module Numeric.GSL.Interpolation , module Data.Complex , setErrorHandlerOff ) where @@ -34,6 +35,7 @@ import Numeric.GSL.Root import Numeric.GSL.ODE import Numeric.GSL.Fitting +import Numeric.GSL.Interpolation import Data.Complex diff --git a/src/Numeric/GSL/Fitting.hs b/src/Numeric/GSL/Fitting.hs --- a/src/Numeric/GSL/Fitting.hs +++ b/src/Numeric/GSL/Fitting.hs @@ -1,3 +1,5 @@ +{-# LANGUAGE FlexibleContexts #-} + {- | Module : Numeric.GSL.Fitting Copyright : (c) Alberto Ruiz 2010 @@ -50,7 +52,7 @@ fitModelScaled, fitModel ) where -import Numeric.LinearAlgebra +import Numeric.LinearAlgebra.HMatrix import Numeric.GSL.Internal import Foreign.Ptr(FunPtr, freeHaskellFunPtr) @@ -80,13 +82,13 @@ nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do - let p = dim xiv - n = dim (f xiv) + let p = size xiv + n = size (f xiv) fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f)) jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac)) rawpath <- createMatrix RowMajor maxit (2+p) - app2 (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) vec xiv mat rawpath "c_nlfit" - let it = round (rawpath @@> (maxit-1,0)) + c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n) # xiv # rawpath #|"c_nlfit" + let it = round (rawpath `atIndex` (maxit-1,0)) path = takeRows it rawpath [sol] = toRows $ dropRows (it-1) path freeHaskellFunPtr fp @@ -99,7 +101,7 @@ ------------------------------------------------------- checkdim1 n _p v - | dim v == n = v + | size v == n = v | otherwise = error $ "Error: "++ show n ++ " components expected in the result of the function supplied to nlFitting" @@ -114,9 +116,9 @@ sol = toList vsol c = max 1 (chi/sqrt (fromIntegral dof)) dof = length dat - (rows cov) - chi = norm2 (fromList $ cost (resMs model) dat sol) + chi = norm_2 (fromList $ cost (resMs model) dat sol) js = fromLists $ jacobian (resDs deriv) dat sol - cov = inv $ trans js <> js + cov = inv $ tr js <> js errs = toList $ scalar c * sqrt (takeDiag cov) diff --git a/src/Numeric/GSL/Fourier.hs b/src/Numeric/GSL/Fourier.hs --- a/src/Numeric/GSL/Fourier.hs +++ b/src/Numeric/GSL/Fourier.hs @@ -1,3 +1,5 @@ +{-# LANGUAGE TypeFamilies #-} + {- | Module : Numeric.GSL.Fourier Copyright : (c) Alberto Ruiz 2006 @@ -16,15 +18,14 @@ ifft ) where -import Data.Packed +import Numeric.LinearAlgebra.HMatrix import Numeric.GSL.Internal -import Data.Complex import Foreign.C.Types import System.IO.Unsafe (unsafePerformIO) genfft code v = unsafePerformIO $ do - r <- createVector (dim v) - app2 (c_fft code) vec v vec r "fft" + r <- createVector (size v) + c_fft code # v # r #|"fft" return r foreign import ccall unsafe "gsl-aux.h fft" c_fft :: CInt -> TCV (TCV Res) @@ -42,3 +43,4 @@ -- | The inverse of 'fft', using /gsl_fft_complex_inverse/. ifft :: Vector (Complex Double) -> Vector (Complex Double) ifft = genfft 1 + diff --git a/src/Numeric/GSL/IO.hs b/src/Numeric/GSL/IO.hs --- a/src/Numeric/GSL/IO.hs +++ b/src/Numeric/GSL/IO.hs @@ -14,7 +14,7 @@ fileDimensions, loadMatrix, fromFile ) where -import Data.Packed +import Numeric.LinearAlgebra.HMatrix hiding(saveMatrix, loadMatrix) import Numeric.GSL.Vector import System.Process(readProcess) diff --git a/src/Numeric/GSL/Internal.hs b/src/Numeric/GSL/Internal.hs --- a/src/Numeric/GSL/Internal.hs +++ b/src/Numeric/GSL/Internal.hs @@ -22,21 +22,20 @@ aux_vTom, createV, createMIO, - module Data.Packed.Development, - check, + module Numeric.LinearAlgebra.Devel, + check,(#),vec, ww2, Res,TV,TM,TCV,TCM ) where -import Data.Packed -import Data.Packed.Development hiding (check) -import Data.Complex +import Numeric.LinearAlgebra.HMatrix +import Numeric.LinearAlgebra.Devel hiding (check) import Foreign.Marshal.Array(copyArray) import Foreign.Ptr(Ptr, FunPtr) import Foreign.C.Types import Foreign.C.String(peekCString) import System.IO.Unsafe(unsafePerformIO) -import Data.Vector.Storable(unsafeWith) +import Data.Vector.Storable as V (unsafeWith,length) import Control.Monad(when) iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double) @@ -87,12 +86,12 @@ createV n fun msg = unsafePerformIO $ do r <- createVector n - app1 fun vec r msg + fun # r #| msg return r createMIO r c fun msg = do res <- createMatrix RowMajor r c - app1 fun mat res msg + fun # res #| msg return res -------------------------------------------------------------------------------- @@ -123,4 +122,16 @@ type TVV = TV (TV Res) type TVM = TV (TM Res) + +ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2 + +vec x f = unsafeWith x $ \p -> do + let v g = do + g (fi $ V.length x) p + f v +{-# INLINE vec #-} + +infixl 1 # +a # b = applyRaw a b +{-# INLINE (#) #-} diff --git a/src/Numeric/GSL/Interpolation.hs b/src/Numeric/GSL/Interpolation.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/GSL/Interpolation.hs @@ -0,0 +1,284 @@ +{- | +Module : Numeric.GSL.Interpolation +Copyright : (c) Matthew Peddie 2015 +License : GPL +Maintainer : Alberto Ruiz +Stability : provisional + +Interpolation routines. + + + +The GSL routines @gsl_spline_eval@ and friends are used, but in spite +of the names, they are not restricted to spline interpolation. The +functions in this module will work for any 'InterpolationMethod'. + +-} + + +module Numeric.GSL.Interpolation ( + -- * Interpolation methods + InterpolationMethod(..) + -- * Evaluation of interpolated functions + , evaluate + , evaluateV + -- * Evaluation of derivatives of interpolated functions + , evaluateDerivative + , evaluateDerivative2 + , evaluateDerivativeV + , evaluateDerivative2V + -- * Evaluation of integrals of interpolated functions + , evaluateIntegral + , evaluateIntegralV +) where + +import Numeric.LinearAlgebra(Vector, fromList, size, Numeric) +import Foreign.C.Types +import Foreign.Marshal.Alloc(alloca) +import Foreign.Ptr(Ptr) +import Foreign.Storable(peek) +import Numeric.GSL.Internal +import System.IO.Unsafe(unsafePerformIO) + +data InterpolationMethod = Linear + | Polynomial + | CSpline + | CSplinePeriodic + | Akima + | AkimaPeriodic + deriving (Eq, Show, Read) + +methodToInt :: Integral a => InterpolationMethod -> a +methodToInt Linear = 0 +methodToInt Polynomial = 1 +methodToInt CSpline = 2 +methodToInt CSplinePeriodic = 3 +methodToInt Akima = 4 +methodToInt AkimaPeriodic = 5 + +dim :: Numeric t => Vector t -> Int +dim = size + +applyCFun hsname cname fun mth xs ys x + | dim xs /= dim ys = error $ + "Error: Vectors of unequal sizes " ++ + show (dim xs) ++ " and " ++ show (dim ys) ++ + " supplied to " ++ hsname + | otherwise = unsafePerformIO $ + flip appVector xs $ \xs' -> + flip appVector ys $ \ys' -> + alloca $ \y' -> do + fun xs' ys' + (fromIntegral $ dim xs) x + (methodToInt mth) y' // check cname + peek y' + +foreign import ccall safe "spline_eval" c_spline_eval + :: Ptr Double -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt + +-------------------------------------------------------------------- +{- | Evaluate a function by interpolating within the given dataset. For +example: + +>>> let xs = vector [1..10] +>>> let ys = vector $ map (**2) [1..10] +>>> evaluateV CSpline xs ys 2.2 +4.818867924528303 + +To successfully @evaluateV xs ys x@, the vectors of corresponding +domain-range values @xs@ and @ys@ must have identical lengths, and +@xs@ must be monotonically increasing. The evaluation point @x@ must +lie between the smallest and largest values in @xs@. + +-} +evaluateV :: InterpolationMethod -- ^ What method to use to interpolate + -> Vector Double -- ^ Data points sampling the domain of the function + -> Vector Double -- ^ Data points sampling the range of the function + -> Double -- ^ Point at which to evaluate the function + -> Double -- ^ Interpolated result +evaluateV = applyCFun "evaluateV" "spline_eval" c_spline_eval + +{- | Evaluate a function by interpolating within the given dataset. For +example: + +>>> let xs = [1..10] +>>> let ys map (**2) [1..10] +>>> evaluate Akima (zip xs ys) 2.2 +4.840000000000001 + +To successfully @evaluate points x@, the domain (@x@) values in +@points@ must be monotonically increasing. The evaluation point @x@ +must lie between the smallest and largest values in the sampled +domain. + +-} +evaluate :: InterpolationMethod -- ^ What method to use to interpolate + -> [(Double, Double)] -- ^ (domain, range) values sampling the function + -> Double -- ^ Point at which to evaluate the function + -> Double -- ^ Interpolated result +evaluate mth pts = + applyCFun "evaluate" "spline_eval" c_spline_eval + mth (fromList xs) (fromList ys) + where + (xs, ys) = unzip pts + +foreign import ccall safe "spline_eval_deriv" c_spline_eval_deriv + :: Ptr Double -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt + +{- | Evaluate the derivative of a function by interpolating within the +given dataset. For example: + +>>> let xs = vector [1..10] +>>> let ys = vector $ map (**2) [1..10] +>>> evaluateDerivativeV CSpline xs ys 2.2 +4.338867924528302 + +To successfully @evaluateDerivativeV xs ys x@, the vectors of +corresponding domain-range values @xs@ and @ys@ must have identical +lengths, and @xs@ must be monotonically increasing. The interpolation +point @x@ must lie between the smallest and largest values in @xs@. + +-} +evaluateDerivativeV :: InterpolationMethod -- ^ What method to use to interpolate + -> Vector Double -- ^ Data points @xs@ sampling the domain of the function + -> Vector Double -- ^ Data points @ys@ sampling the range of the function + -> Double -- ^ Point @x@ at which to evaluate the derivative + -> Double -- ^ Interpolated result +evaluateDerivativeV = + applyCFun "evaluateDerivativeV" "spline_eval_deriv" c_spline_eval_deriv + +{- | Evaluate the derivative of a function by interpolating within the +given dataset. For example: + +>>> let xs = [1..10] +>>> let ys map (**2) [1..10] +>>> evaluateDerivative Akima (zip xs ys) 2.2 +4.4 + +To successfully @evaluateDerivative points x@, the domain (@x@) values +in @points@ must be monotonically increasing. The evaluation point +@x@ must lie between the smallest and largest values in the sampled +domain. + +-} +evaluateDerivative :: InterpolationMethod -- ^ What method to use to interpolate + -> [(Double, Double)] -- ^ (domain, range) points sampling the function + -> Double -- ^ Point @x@ at which to evaluate the derivative + -> Double -- ^ Interpolated result +evaluateDerivative mth pts = + applyCFun "evaluateDerivative" "spline_eval_deriv" c_spline_eval_deriv + mth (fromList xs) (fromList ys) + where + (xs, ys) = unzip pts + +foreign import ccall safe "spline_eval_deriv2" c_spline_eval_deriv2 + :: Ptr Double -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt + +{- | Evaluate the second derivative of a function by interpolating within the +given dataset. For example: + +>>> let xs = vector [1..10] +>>> let ys = vector $ map (**2) [1..10] +>>> evaluateDerivative2V CSpline xs ys 2.2 +2.4 + +To successfully @evaluateDerivative2V xs ys x@, the vectors @xs@ and +@ys@ must have identical lengths, and @xs@ must be monotonically +increasing. The evaluation point @x@ must lie between the smallest +and largest values in @xs@. + +-} +evaluateDerivative2V :: InterpolationMethod -- ^ What method to use to interpolate + -> Vector Double -- ^ Data points @xs@ sampling the domain of the function + -> Vector Double -- ^ Data points @ys@ sampling the range of the function + -> Double -- ^ Point @x@ at which to evaluate the second derivative + -> Double -- ^ Interpolated result +evaluateDerivative2V = + applyCFun "evaluateDerivative2V" "spline_eval_deriv2" c_spline_eval_deriv2 + +{- | Evaluate the second derivative of a function by interpolating +within the given dataset. For example: + +>>> let xs = [1..10] +>>> let ys map (**2) [1..10] +>>> evaluateDerivative2 Akima (zip xs ys) 2.2 +2.0 + +To successfully @evaluateDerivative2 points x@, the domain (@x@) +values in @points@ must be monotonically increasing. The evaluation +point @x@ must lie between the smallest and largest values in the +sampled domain. + +-} +evaluateDerivative2 :: InterpolationMethod -- ^ What method to use to interpolate + -> [(Double, Double)] -- ^ (domain, range) points sampling the function + -> Double -- ^ Point @x@ at which to evaluate the second derivative + -> Double -- ^ Interpolated result +evaluateDerivative2 mth pts = + applyCFun "evaluateDerivative2" "spline_eval_deriv2" c_spline_eval_deriv2 + mth (fromList xs) (fromList ys) + where + (xs, ys) = unzip pts + +foreign import ccall safe "spline_eval_integ" c_spline_eval_integ + :: Ptr Double -> Ptr Double -> CUInt -> Double -> Double -> CInt -> Ptr Double -> IO CInt + +applyCIntFun hsname cname fun mth xs ys a b + | dim xs /= dim ys = error $ + "Error: Vectors of unequal sizes " ++ + show (dim xs) ++ " and " ++ show (dim ys) ++ + " supplied to " ++ hsname + | otherwise = unsafePerformIO $ + flip appVector xs $ \xs' -> + flip appVector ys $ \ys' -> + alloca $ \y' -> do + fun xs' ys' + (fromIntegral $ dim xs) a b + (methodToInt mth) y' // check cname + peek y' + +{- | Evaluate the definite integral of a function by interpolating +within the given dataset. For example: + +>>> let xs = vector [1..10] +>>> let ys = vector $ map (**2) [1..10] +>>> evaluateIntegralV CSpline xs ys 2.2 5.5 +51.89853207547169 + +To successfully @evaluateIntegralV xs ys a b@, the vectors @xs@ and +@ys@ must have identical lengths, and @xs@ must be monotonically +increasing. The integration bounds @a@ and @b@ must lie between the +smallest and largest values in @xs@. + +-} +evaluateIntegralV :: InterpolationMethod -- ^ What method to use to interpolate + -> Vector Double -- ^ Data points @xs@ sampling the domain of the function + -> Vector Double -- ^ Data points @ys@ sampling the range of the function + -> Double -- ^ Lower integration bound @a@ + -> Double -- ^ Upper integration bound @b@ + -> Double -- ^ Resulting area +evaluateIntegralV = + applyCIntFun "evaluateIntegralV" "spline_eval_integ" c_spline_eval_integ + +{- | Evaluate the definite integral of a function by interpolating +within the given dataset. For example: + +>>> let xs = [1..10] +>>> let ys = map (**2) [1..10] +>>> evaluateIntegralV CSpline (zip xs ys) (2.2, 5.5) +51.909 + +To successfully @evaluateIntegral points (a, b)@, the domain (@x@) +values of @points@ must be monotonically increasing. The integration +bounds @a@ and @b@ must lie between the smallest and largest values in +the sampled domain.. +-} +evaluateIntegral :: InterpolationMethod -- ^ What method to use to interpolate + -> [(Double, Double)] -- ^ (domain, range) points sampling the function + -> (Double, Double) -- ^ Integration bounds (@a@, @b@) + -> Double -- ^ Resulting area +evaluateIntegral mth pts (a, b) = + applyCIntFun "evaluateIntegral" "spline_eval_integ" c_spline_eval_integ + mth (fromList xs) (fromList ys) a b + where + (xs, ys) = unzip pts diff --git a/src/Numeric/GSL/LinearAlgebra.hs b/src/Numeric/GSL/LinearAlgebra.hs --- a/src/Numeric/GSL/LinearAlgebra.hs +++ b/src/Numeric/GSL/LinearAlgebra.hs @@ -15,7 +15,7 @@ fileDimensions, loadMatrix, fromFile ) where -import Data.Packed +import Numeric.LinearAlgebra.HMatrix hiding (RandDist,randomVector,saveMatrix,loadMatrix) import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) import Foreign.Marshal.Alloc(free) @@ -40,7 +40,7 @@ -> Vector Double randomVector seed dist n = unsafePerformIO $ do r <- createVector n - app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector" + c_random_vector (fi seed) ((fi.fromEnum) dist) # r #|"randomVector" return r foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV @@ -56,7 +56,7 @@ charname <- newCString filename charfmt <- newCString fmt let o = if orderOf m == RowMajor then 1 else 0 - app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf" + matrix_fprintf charname charfmt o # m #|"matrix_fprintf" free charname free charfmt @@ -69,7 +69,7 @@ fscanfVector filename n = do charname <- newCString filename res <- createVector n - app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf" + gsl_vector_fscanf charname # res #|"gsl_vector_fscanf" free charname return res @@ -80,7 +80,7 @@ fprintfVector filename fmt v = do charname <- newCString filename charfmt <- newCString fmt - app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf" + gsl_vector_fprintf charname charfmt # v #|"gsl_vector_fprintf" free charname free charfmt @@ -91,7 +91,7 @@ freadVector filename n = do charname <- newCString filename res <- createVector n - app1 (gsl_vector_fread charname) vec res "gsl_vector_fread" + gsl_vector_fread charname # res #| "gsl_vector_fread" free charname return res @@ -101,7 +101,7 @@ fwriteVector :: FilePath -> Vector Double -> IO () fwriteVector filename v = do charname <- newCString filename - app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite" + gsl_vector_fwrite charname # v #|"gsl_vector_fwrite" free charname foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV diff --git a/src/Numeric/GSL/Minimization.hs b/src/Numeric/GSL/Minimization.hs --- a/src/Numeric/GSL/Minimization.hs +++ b/src/Numeric/GSL/Minimization.hs @@ -1,3 +1,6 @@ +{-# LANGUAGE FlexibleContexts #-} + + {- | Module : Numeric.GSL.Minimization Copyright : (c) Alberto Ruiz 2006-9 @@ -56,7 +59,7 @@ ) where -import Data.Packed +import Numeric.LinearAlgebra.HMatrix hiding(step) import Numeric.GSL.Internal import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr) @@ -99,7 +102,7 @@ rawpath <- createMIO maxit 4 (c_uniMinize m fp epsrel (fi maxit) xmin xl xu) "uniMinimize" - let it = round (rawpath @@> (maxit-1,0)) + let it = round (rawpath `atIndex` (maxit-1,0)) path = takeRows it rawpath [sol] = toLists $ dropRows (it-1) path freeHaskellFunPtr fp @@ -134,16 +137,16 @@ minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList sz) (f.toList) (fromList xi) where v2l (v,m) = (toList v, m) -ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2 + minimizeV method eps maxit szv f xiv = unsafePerformIO $ do - let n = dim xiv + let n = size xiv fp <- mkVecfun (iv f) rawpath <- ww2 vec xiv vec szv $ \xiv' szv' -> createMIO maxit (n+3) (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv') "minimize" - let it = round (rawpath @@> (maxit-1,0)) + let it = round (rawpath `atIndex` (maxit-1,0)) path = takeRows it rawpath sol = flatten $ dropColumns 3 $ dropRows (it-1) path freeHaskellFunPtr fp @@ -191,7 +194,7 @@ minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do - let n = dim xiv + let n = size xiv f' = f df' = (checkdim1 n . df) fp <- mkVecfun (iv f') @@ -200,7 +203,7 @@ createMIO maxit (n+2) (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv') "minimizeD" - let it = round (rawpath @@> (maxit-1,0)) + let it = round (rawpath `atIndex` (maxit-1,0)) path = takeRows it rawpath sol = flatten $ dropColumns 2 $ dropRows (it-1) path freeHaskellFunPtr fp @@ -217,6 +220,6 @@ --------------------------------------------------------------------- checkdim1 n v - | dim v == n = v + | size v == n = v | otherwise = error $ "Error: "++ show n ++ " components expected in the result of the gradient supplied to minimizeD" diff --git a/src/Numeric/GSL/ODE.hs b/src/Numeric/GSL/ODE.hs --- a/src/Numeric/GSL/ODE.hs +++ b/src/Numeric/GSL/ODE.hs @@ -1,3 +1,6 @@ +{-# LANGUAGE FlexibleContexts #-} + + {- | Module : Numeric.GSL.ODE Copyright : (c) Alberto Ruiz 2010 @@ -29,10 +32,10 @@ ----------------------------------------------------------------------------- module Numeric.GSL.ODE ( - odeSolve, odeSolveV, ODEMethod(..), Jacobian + odeSolve, odeSolveV, odeSolveVWith, ODEMethod(..), Jacobian, StepControl(..) ) where -import Data.Packed +import Numeric.LinearAlgebra.HMatrix import Numeric.GSL.Internal import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr) @@ -41,9 +44,10 @@ ------------------------------------------------------------------------- -type TVV = TV (TV Res) -type TVM = TV (TM Res) -type TVVM = TV (TV (TM Res)) +type TVV = TV (TV Res) +type TVM = TV (TM Res) +type TVVM = TV (TV (TM Res)) +type TVVVM = TV (TV (TV (TM Res))) type Jacobian = Double -> Vector Double -> Matrix Double @@ -60,73 +64,105 @@ | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12. | MSBDF Jacobian -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems. +-- | Adaptive step-size control functions +data StepControl = X Double Double -- ^ abs. and rel. tolerance for x(t) + | X' Double Double -- ^ abs. and rel. tolerance for x'(t) + | XX' Double Double Double Double -- ^ include both via rel. tolerance scaling factors a_x, a_x' + | ScXX' Double Double Double Double (Vector Double) -- ^ scale abs. tolerance of x(t) components -- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists. odeSolve - :: (Double -> [Double] -> [Double]) -- ^ xdot(t,x) + :: (Double -> [Double] -> [Double]) -- ^ x'(t,x) -> [Double] -- ^ initial conditions -> Vector Double -- ^ desired solution times -> Matrix Double -- ^ solution odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts - where hi = (ts@>1 - ts@>0)/100 + where hi = (ts!1 - ts!0)/100 epsAbs = 1.49012e-08 - epsRel = 1.49012e-08 - l2v f = \t -> fromList . f t . toList + epsRel = epsAbs + l2v f = \t -> fromList . f t . toList --- | Evolution of the system with adaptive step-size control. +-- | A version of 'odeSolveVWith' with reasonable default step control. odeSolveV :: ODEMethod - -> Double -- ^ initial step size - -> Double -- ^ absolute tolerance for the state vector - -> Double -- ^ relative tolerance for the state vector - -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) + -> Double -- ^ initial step size + -> Double -- ^ absolute tolerance for the state vector + -> Double -- ^ relative tolerance for the state vector + -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x) -> Vector Double -- ^ initial conditions -> Vector Double -- ^ desired solution times -> Matrix Double -- ^ solution -odeSolveV RK2 = odeSolveV' 0 Nothing -odeSolveV RK4 = odeSolveV' 1 Nothing -odeSolveV RKf45 = odeSolveV' 2 Nothing -odeSolveV RKck = odeSolveV' 3 Nothing -odeSolveV RK8pd = odeSolveV' 4 Nothing -odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac) -odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac) -odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac) -odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac) -odeSolveV MSAdams = odeSolveV' 9 Nothing -odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac) - +odeSolveV meth hi epsAbs epsRel = odeSolveVWith meth (XX' epsAbs epsRel 1 1) hi -odeSolveV' - :: CInt - -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian - -> Double -- ^ initial step size - -> Double -- ^ absolute tolerance for the state vector - -> Double -- ^ relative tolerance for the state vector - -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x) +-- | Evolution of the system with adaptive step-size control. +odeSolveVWith + :: ODEMethod + -> StepControl + -> Double -- ^ initial step size + -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x) -> Vector Double -- ^ initial conditions -> Vector Double -- ^ desired solution times -> Matrix Double -- ^ solution -odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do - let n = dim xiv - fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) - jp <- case mbjac of - Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t)) - Nothing -> return nullFunPtr - sol <- vec xiv $ \xiv' -> - vec (checkTimes ts) $ \ts' -> - createMIO (dim ts) n - (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' ) - "ode" - freeHaskellFunPtr fp - return sol +odeSolveVWith method control = odeSolveVWith' m mbj c epsAbs epsRel aX aX' mbsc + where (m, mbj) = case method of + RK2 -> (0 , Nothing ) + RK4 -> (1 , Nothing ) + RKf45 -> (2 , Nothing ) + RKck -> (3 , Nothing ) + RK8pd -> (4 , Nothing ) + RK2imp jac -> (5 , Just jac) + RK4imp jac -> (6 , Just jac) + BSimp jac -> (7 , Just jac) + RK1imp jac -> (8 , Just jac) + MSAdams -> (9 , Nothing ) + MSBDF jac -> (10, Just jac) + (c, epsAbs, epsRel, aX, aX', mbsc) = case control of + X ea er -> (0, ea, er, 1 , 0 , Nothing) + X' ea er -> (0, ea, er, 0 , 1 , Nothing) + XX' ea er ax ax' -> (0, ea, er, ax, ax', Nothing) + ScXX' ea er ax ax' sc -> (1, ea, er, ax, ax', Just sc) +odeSolveVWith' + :: CInt -- ^ stepping function + -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian + -> CInt -- ^ step-size control function + -> Double -- ^ absolute tolerance for step-size control + -> Double -- ^ relative tolerance for step-size control + -> Double -- ^ scaling factor for relative tolerance of x(t) + -> Double -- ^ scaling factor for relative tolerance of x'(t) + -> Maybe (Vector Double) -- ^ optional scaling for absolute error + -> Double -- ^ initial step size + -> (Double -> Vector Double -> Vector Double) -- ^ x'(t,x) + -> Vector Double -- ^ initial conditions + -> Vector Double -- ^ desired solution times + -> Matrix Double -- ^ solution +odeSolveVWith' method mbjac control epsAbs epsRel aX aX' mbsc h f xiv ts = + unsafePerformIO $ do + let n = size xiv + sc = case mbsc of + Just scv -> checkdim1 n scv + Nothing -> xiv + fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) + jp <- case mbjac of + Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t)) + Nothing -> return nullFunPtr + sol <- vec sc $ \sc' -> vec xiv $ \xiv' -> + vec (checkTimes ts) $ \ts' -> createMIO (size ts) n + (ode_c method control h epsAbs epsRel aX aX' fp jp + // sc' // xiv' // ts' ) + "ode" + freeHaskellFunPtr fp + return sol + foreign import ccall safe "ode" - ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM + ode_c :: CInt -> CInt -> Double + -> Double -> Double -> Double -> Double + -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVVM ------------------------------------------------------- checkdim1 n v - | dim v == n = v + | size v == n = v | otherwise = error $ "Error: "++ show n ++ " components expected in the result of the function supplied to odeSolve" @@ -135,6 +171,6 @@ | otherwise = error $ "Error: "++ show n ++ "x" ++ show n ++ " Jacobian expected in odeSolve" -checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts +checkTimes ts | size ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts | otherwise = error "odeSolve requires increasing times" where ts' = toList ts diff --git a/src/Numeric/GSL/Polynomials.hs b/src/Numeric/GSL/Polynomials.hs --- a/src/Numeric/GSL/Polynomials.hs +++ b/src/Numeric/GSL/Polynomials.hs @@ -16,9 +16,8 @@ polySolve ) where -import Data.Packed +import Numeric.LinearAlgebra.HMatrix import Numeric.GSL.Internal -import Data.Complex import System.IO.Unsafe (unsafePerformIO) #if __GLASGOW_HASKELL__ >= 704 @@ -47,9 +46,9 @@ polySolve = toList . polySolve' . fromList polySolve' :: Vector Double -> Vector (Complex Double) -polySolve' v | dim v > 1 = unsafePerformIO $ do - r <- createVector (dim v-1) - app2 c_polySolve vec v vec r "polySolve" +polySolve' v | size v > 1 = unsafePerformIO $ do + r <- createVector (size v-1) + c_polySolve # v # r #| "polySolve" return r | otherwise = error "polySolve on a polynomial of degree zero" diff --git a/src/Numeric/GSL/Random.hs b/src/Numeric/GSL/Random.hs --- a/src/Numeric/GSL/Random.hs +++ b/src/Numeric/GSL/Random.hs @@ -21,11 +21,13 @@ ) where import Numeric.GSL.Vector -import Numeric.LinearAlgebra(cholSH) -import Numeric.Container hiding ( +import Numeric.LinearAlgebra.HMatrix hiding ( randomVector, gaussianSample, - uniformSample + uniformSample, + Seed, + rand, + randn ) import System.Random(randomIO) @@ -40,10 +42,10 @@ -> Matrix Double -- ^ covariance matrix -> Matrix Double -- ^ result gaussianSample seed n med cov = m where - c = dim med + c = size med meds = konst 1 n `outer` med rs = reshape c $ randomVector seed Gaussian (c * n) - m = rs `mXm` cholSH cov `add` meds + m = rs <> cholSH cov + meds -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate -- uniform distribution. @@ -55,10 +57,10 @@ (as,bs) = unzip rgs a = fromList as cs = zipWith subtract as bs - d = dim a + d = size a dat = toRows $ reshape n $ randomVector seed Uniform (n*d) am = konst 1 n `outer` a - m = fromColumns (zipWith scale cs dat) `add` am + m = fromColumns (zipWith scale cs dat) + am -- | pseudorandom matrix with uniform elements between 0 and 1 randm :: RandDist diff --git a/src/Numeric/GSL/Root.hs b/src/Numeric/GSL/Root.hs --- a/src/Numeric/GSL/Root.hs +++ b/src/Numeric/GSL/Root.hs @@ -1,3 +1,5 @@ +{-# LANGUAGE FlexibleContexts #-} + {- | Module : Numeric.GSL.Root Copyright : (c) Alberto Ruiz 2009 @@ -39,7 +41,7 @@ rootJ, RootMethodJ(..), ) where -import Data.Packed +import Numeric.LinearAlgebra.HMatrix import Numeric.GSL.Internal import Foreign.Ptr(FunPtr, freeHaskellFunPtr) import Foreign.C.Types @@ -69,7 +71,7 @@ rawpath <- createMIO maxit 4 (c_root m fp epsrel (fi maxit) xl xu) "root" - let it = round (rawpath @@> (maxit-1,0)) + let it = round (rawpath `atIndex` (maxit-1,0)) path = takeRows it rawpath [sol] = toLists $ dropRows (it-1) path freeHaskellFunPtr fp @@ -100,7 +102,7 @@ rawpath <- createMIO maxit 2 (c_rootj m fp dfp epsrel (fi maxit) x) "rootj" - let it = round (rawpath @@> (maxit-1,0)) + let it = round (rawpath `atIndex` (maxit-1,0)) path = takeRows it rawpath [sol] = toLists $ dropRows (it-1) path freeHaskellFunPtr fp @@ -132,13 +134,13 @@ rootGen m f xi epsabs maxit = unsafePerformIO $ do let xiv = fromList xi - n = dim xiv + n = size xiv fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) rawpath <- vec xiv $ \xiv' -> createMIO maxit (2*n+1) (c_multiroot m fp epsabs (fi maxit) // xiv') "multiroot" - let it = round (rawpath @@> (maxit-1,0)) + let it = round (rawpath `atIndex` (maxit-1,0)) path = takeRows it rawpath [sol] = toLists $ dropRows (it-1) path freeHaskellFunPtr fp @@ -169,14 +171,14 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do let xiv = fromList xi - n = dim xiv + n = size xiv fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList)) rawpath <- vec xiv $ \xiv' -> createMIO maxit (2*n+1) (c_multirootj m fp jp epsabs (fi maxit) // xiv') "multiroot" - let it = round (rawpath @@> (maxit-1,0)) + let it = round (rawpath `atIndex` (maxit-1,0)) path = takeRows it rawpath [sol] = toLists $ dropRows (it-1) path freeHaskellFunPtr fp @@ -189,7 +191,7 @@ ------------------------------------------------------- checkdim1 n v - | dim v == n = v + | size v == n = v | otherwise = error $ "Error: "++ show n ++ " components expected in the result of the function supplied to root" diff --git a/src/Numeric/GSL/SimulatedAnnealing.hs b/src/Numeric/GSL/SimulatedAnnealing.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/GSL/SimulatedAnnealing.hs @@ -0,0 +1,245 @@ +{- | +Module : Numeric.GSL.Interpolation +Copyright : (c) Matthew Peddie 2015 +License : GPL +Maintainer : Alberto Ruiz +Stability : provisional + +Simulated annealing routines. + + + +Here is a translation of the simple example given in +: + +> import Numeric.GSL.SimulatedAnnealing +> import Numeric.LinearAlgebra.HMatrix +> +> main = print $ simanSolve 0 1 exampleParams 15.5 exampleE exampleM exampleS (Just show) +> +> exampleParams = SimulatedAnnealingParams 200 1000 1.0 1.0 0.008 1.003 2.0e-6 +> +> exampleE x = exp (-(x - 1)**2) * sin (8 * x) +> +> exampleM x y = abs $ x - y +> +> exampleS rands stepSize current = (rands ! 0) * 2 * stepSize - stepSize + current + +The manual states: + +> The first example, in one dimensional Cartesian space, sets up an +> energy function which is a damped sine wave; this has many local +> minima, but only one global minimum, somewhere between 1.0 and +> 1.5. The initial guess given is 15.5, which is several local minima +> away from the global minimum. + +This global minimum is around 1.36. + +-} +{-# OPTIONS_GHC -Wall #-} + +module Numeric.GSL.SimulatedAnnealing ( + -- * Searching for minima + simanSolve + -- * Configuring the annealing process + , SimulatedAnnealingParams(..) + ) where + +import Numeric.GSL.Internal +import Numeric.LinearAlgebra.HMatrix hiding(step) + +import Data.Vector.Storable(generateM) +import Foreign.Storable(Storable(..)) +import Foreign.Marshal.Utils(with) +import Foreign.Ptr(Ptr, FunPtr, nullFunPtr) +import Foreign.StablePtr(StablePtr, newStablePtr, deRefStablePtr, freeStablePtr) +import Foreign.C.Types +import System.IO.Unsafe(unsafePerformIO) + +import System.IO (hFlush, stdout) + +import Data.IORef (IORef, newIORef, writeIORef, readIORef, modifyIORef') + +-- | 'SimulatedAnnealingParams' is a translation of the +-- @gsl_siman_params_t@ structure documented in +-- , +-- which controls the simulated annealing algorithm. +-- +-- The annealing process is parameterized by the Boltzmann +-- distribution and the /cooling schedule/. For more details, see +-- . +data SimulatedAnnealingParams = SimulatedAnnealingParams { + n_tries :: CInt -- ^ The number of points to try for each step. + , iters_fixed_T :: CInt -- ^ The number of iterations at each temperature + , step_size :: Double -- ^ The maximum step size in the random walk + , boltzmann_k :: Double -- ^ Boltzmann distribution parameter + , cooling_t_initial :: Double -- ^ Initial temperature + , cooling_mu_t :: Double -- ^ Cooling rate parameter + , cooling_t_min :: Double -- ^ Final temperature + } deriving (Eq, Show, Read) + +instance Storable SimulatedAnnealingParams where + sizeOf p = sizeOf (n_tries p) + + sizeOf (iters_fixed_T p) + + sizeOf (step_size p) + + sizeOf (boltzmann_k p) + + sizeOf (cooling_t_initial p) + + sizeOf (cooling_mu_t p) + + sizeOf (cooling_t_min p) + -- TODO(MP): is this safe? + alignment p = alignment (step_size p) + -- TODO(MP): Is there a more automatic way to write these? + peek ptr = SimulatedAnnealingParams <$> + peekByteOff ptr 0 <*> + peekByteOff ptr i <*> + peekByteOff ptr (2*i) <*> + peekByteOff ptr (2*i + d) <*> + peekByteOff ptr (2*i + 2*d) <*> + peekByteOff ptr (2*i + 3*d) <*> + peekByteOff ptr (2*i + 4*d) + where + i = sizeOf (0 :: CInt) + d = sizeOf (0 :: Double) + poke ptr sap = do + pokeByteOff ptr 0 (n_tries sap) + pokeByteOff ptr i (iters_fixed_T sap) + pokeByteOff ptr (2*i) (step_size sap) + pokeByteOff ptr (2*i + d) (boltzmann_k sap) + pokeByteOff ptr (2*i + 2*d) (cooling_t_initial sap) + pokeByteOff ptr (2*i + 3*d) (cooling_mu_t sap) + pokeByteOff ptr (2*i + 4*d) (cooling_t_min sap) + where + i = sizeOf (0 :: CInt) + d = sizeOf (0 :: Double) + +-- We use a StablePtr to an IORef so that we can keep hold of +-- StablePtr values but mutate their contents. A simple 'StablePtr a' +-- won't work, since we'd have no way to write 'copyConfig'. +type P a = StablePtr (IORef a) + +copyConfig :: P a -> P a -> IO () +copyConfig src' dest' = do + dest <- deRefStablePtr dest' + src <- deRefStablePtr src' + readIORef src >>= writeIORef dest + +copyConstructConfig :: P a -> IO (P a) +copyConstructConfig x = do + conf <- deRefRead x + newconf <- newIORef conf + newStablePtr newconf + +destroyConfig :: P a -> IO () +destroyConfig p = do + freeStablePtr p + +deRefRead :: P a -> IO a +deRefRead p = deRefStablePtr p >>= readIORef + +wrapEnergy :: (a -> Double) -> P a -> Double +wrapEnergy f p = unsafePerformIO $ f <$> deRefRead p + +wrapMetric :: (a -> a -> Double) -> P a -> P a -> Double +wrapMetric f x y = unsafePerformIO $ f <$> deRefRead x <*> deRefRead y + +wrapStep :: Int + -> (Vector Double -> Double -> a -> a) + -> GSLRNG + -> P a + -> Double + -> IO () +wrapStep nrand f (GSLRNG rng) confptr stepSize = do + v <- generateM nrand (\_ -> gslRngUniform rng) + conf <- deRefStablePtr confptr + modifyIORef' conf $ f v stepSize + +wrapPrint :: (a -> String) -> P a -> IO () +wrapPrint pf ptr = deRefRead ptr >>= putStr . pf >> hFlush stdout + +foreign import ccall safe "wrapper" + mkEnergyFun :: (P a -> Double) -> IO (FunPtr (P a -> Double)) + +foreign import ccall safe "wrapper" + mkMetricFun :: (P a -> P a -> Double) -> IO (FunPtr (P a -> P a -> Double)) + +foreign import ccall safe "wrapper" + mkStepFun :: (GSLRNG -> P a -> Double -> IO ()) + -> IO (FunPtr (GSLRNG -> P a -> Double -> IO ())) + +foreign import ccall safe "wrapper" + mkCopyFun :: (P a -> P a -> IO ()) -> IO (FunPtr (P a -> P a -> IO ())) + +foreign import ccall safe "wrapper" + mkCopyConstructorFun :: (P a -> IO (P a)) -> IO (FunPtr (P a -> IO (P a))) + +foreign import ccall safe "wrapper" + mkDestructFun :: (P a -> IO ()) -> IO (FunPtr (P a -> IO ())) + +newtype GSLRNG = GSLRNG (Ptr GSLRNG) + +foreign import ccall safe "gsl_rng.h gsl_rng_uniform" + gslRngUniform :: Ptr GSLRNG -> IO Double + +foreign import ccall safe "gsl-aux.h siman" + siman :: CInt -- ^ RNG seed (for repeatability) + -> Ptr SimulatedAnnealingParams -- ^ params + -> P a -- ^ Configuration + -> FunPtr (P a -> Double) -- ^ Energy functional + -> FunPtr (P a -> P a -> Double) -- ^ Metric definition + -> FunPtr (GSLRNG -> P a -> Double -> IO ()) -- ^ Step evaluation + -> FunPtr (P a -> P a -> IO ()) -- ^ Copy config + -> FunPtr (P a -> IO (P a)) -- ^ Copy constructor for config + -> FunPtr (P a -> IO ()) -- ^ Destructor for config + -> FunPtr (P a -> IO ()) -- ^ Print function + -> IO CInt + +-- | +-- Calling +-- +-- > simanSolve seed nrand params x0 e m step print +-- +-- performs a simulated annealing search through a given space. So +-- that any configuration type may be used, the space is specified by +-- providing the functions @e@ (the energy functional) and @m@ (the +-- metric definition). @x0@ is the initial configuration of the +-- system. The simulated annealing steps are generated using the +-- user-provided function @step@, which should randomly construct a +-- new system configuration. +-- +-- If 'Nothing' is passed instead of a printing function, no +-- incremental output will be generated. Otherwise, the GSL-formatted +-- output, including the configuration description the user function +-- generates, will be printed to stdout. +-- +-- Each time the step function is called, it is supplied with a random +-- vector containing @nrand@ 'Double' values, uniformly distributed in +-- @[0, 1)@. It should use these values to generate its new +-- configuration. +simanSolve :: Int -- ^ Seed for the random number generator + -> Int -- ^ @nrand@, the number of random 'Double's the + -- step function requires + -> SimulatedAnnealingParams -- ^ Parameters to configure the solver + -> a -- ^ Initial configuration @x0@ + -> (a -> Double) -- ^ Energy functional @e@ + -> (a -> a -> Double) -- ^ Metric definition @m@ + -> (Vector Double -> Double -> a -> a) -- ^ Stepping function @step@ + -> Maybe (a -> String) -- ^ Optional printing function + -> a -- ^ Best configuration the solver has found +simanSolve seed nrand params conf e m step printfun = + unsafePerformIO $ with params $ \paramptr -> do + ewrap <- mkEnergyFun $ wrapEnergy e + mwrap <- mkMetricFun $ wrapMetric m + stepwrap <- mkStepFun $ wrapStep nrand step + confptr <- newIORef conf >>= newStablePtr + cpwrap <- mkCopyFun copyConfig + ccwrap <- mkCopyConstructorFun copyConstructConfig + dwrap <- mkDestructFun destroyConfig + pwrap <- case printfun of + Nothing -> return nullFunPtr + Just pf -> mkDestructFun $ wrapPrint pf + siman (fromIntegral seed) + paramptr confptr + ewrap mwrap stepwrap cpwrap ccwrap dwrap pwrap // check "siman" + result <- deRefRead confptr + freeStablePtr confptr + return result diff --git a/src/Numeric/GSL/Vector.hs b/src/Numeric/GSL/Vector.hs --- a/src/Numeric/GSL/Vector.hs +++ b/src/Numeric/GSL/Vector.hs @@ -14,8 +14,7 @@ fwriteVector, freadVector, fprintfVector, fscanfVector ) where -import Data.Packed -import Numeric.LinearAlgebra(RandDist(..)) +import Numeric.LinearAlgebra.HMatrix hiding(randomVector, saveMatrix) import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) import Foreign.Marshal.Alloc(free) @@ -35,7 +34,7 @@ -> Vector Double randomVector seed dist n = unsafePerformIO $ do r <- createVector n - app1 (c_random_vector_GSL (fi seed) ((fi.fromEnum) dist)) vec r "randomVectorGSL" + c_random_vector_GSL (fi seed) ((fi.fromEnum) dist) # r #|"randomVectorGSL" return r foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: CInt -> CInt -> TV @@ -51,7 +50,7 @@ charname <- newCString filename charfmt <- newCString fmt let o = if orderOf m == RowMajor then 1 else 0 - app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf" + matrix_fprintf charname charfmt o # m #|"matrix_fprintf" free charname free charfmt @@ -64,7 +63,7 @@ fscanfVector filename n = do charname <- newCString filename res <- createVector n - app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf" + gsl_vector_fscanf charname # res #|"gsl_vector_fscanf" free charname return res @@ -75,7 +74,7 @@ fprintfVector filename fmt v = do charname <- newCString filename charfmt <- newCString fmt - app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf" + gsl_vector_fprintf charname charfmt # v #|"gsl_vector_fprintf" free charname free charfmt @@ -86,7 +85,7 @@ freadVector filename n = do charname <- newCString filename res <- createVector n - app1 (gsl_vector_fread charname) vec res "gsl_vector_fread" + gsl_vector_fread charname # res #|"gsl_vector_fread" free charname return res @@ -96,7 +95,7 @@ fwriteVector :: FilePath -> Vector Double -> IO () fwriteVector filename v = do charname <- newCString filename - app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite" + gsl_vector_fwrite charname # v #|"gsl_vector_fwrite" free charname foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV diff --git a/src/Numeric/GSL/gsl-aux.c b/src/Numeric/GSL/gsl-aux.c --- a/src/Numeric/GSL/gsl-aux.c +++ b/src/Numeric/GSL/gsl-aux.c @@ -34,7 +34,10 @@ #include #include #include +#include #include +#include + #include #include @@ -140,7 +143,118 @@ return 0; } +int spline_eval(const double xa[], const double ya[], unsigned int size, + double x, int method, double *y) { + DEBUGMSG("spline_eval"); + const gsl_interp_type *T; + switch (method) { + case 0: { T = gsl_interp_linear; break; } + case 1: { T = gsl_interp_polynomial; break; } + case 2: { T = gsl_interp_cspline; break; } + case 3: { T = gsl_interp_cspline_periodic; break; } + case 4: { T = gsl_interp_akima; break; } + case 5: { T = gsl_interp_akima_periodic; break; } + default: ERROR(BAD_CODE); + } + gsl_spline *spline = gsl_spline_alloc(T, size); + if (NULL == spline) ERROR(MEM); + const int initres = gsl_spline_init(spline, xa, ya, size); + CHECK(initres,initres); + gsl_interp_accel *acc = gsl_interp_accel_alloc(); + if (NULL == acc) { gsl_spline_free(spline); ERROR(MEM); }; + + const int res = gsl_spline_eval_e(spline, x, acc, y); + CHECK(res,res); + gsl_interp_accel_free(acc); + gsl_spline_free(spline); + OK +} + +int spline_eval_deriv(const double xa[], const double ya[], unsigned int size, + double x, int method, double *y) { + DEBUGMSG("spline_eval_deriv"); + const gsl_interp_type *T; + switch (method) { + case 0: { T = gsl_interp_linear; break; } + case 1: { T = gsl_interp_polynomial; break; } + case 2: { T = gsl_interp_cspline; break; } + case 3: { T = gsl_interp_cspline_periodic; break; } + case 4: { T = gsl_interp_akima; break; } + case 5: { T = gsl_interp_akima_periodic; break; } + default: ERROR(BAD_CODE); + } + + gsl_spline *spline = gsl_spline_alloc(T, size); + if (NULL == spline) ERROR(MEM); + const int initres = gsl_spline_init(spline, xa, ya, size); + CHECK(initres,initres); + gsl_interp_accel *acc = gsl_interp_accel_alloc(); + if (NULL == acc) { gsl_spline_free(spline); ERROR(MEM); }; + + const int res = gsl_spline_eval_deriv_e(spline, x, acc, y); + CHECK(res,res); + gsl_interp_accel_free(acc); + gsl_spline_free(spline); + OK +} + +int spline_eval_deriv2(const double xa[], const double ya[], unsigned int size, + double x, int method, double *y) { + DEBUGMSG("spline_eval_deriv2"); + const gsl_interp_type *T; + switch (method) { + case 0: { T = gsl_interp_linear; break; } + case 1: { T = gsl_interp_polynomial; break; } + case 2: { T = gsl_interp_cspline; break; } + case 3: { T = gsl_interp_cspline_periodic; break; } + case 4: { T = gsl_interp_akima; break; } + case 5: { T = gsl_interp_akima_periodic; break; } + default: ERROR(BAD_CODE); + } + + gsl_spline *spline = gsl_spline_alloc(T, size); + if (NULL == spline) ERROR(MEM); + const int initres = gsl_spline_init(spline, xa, ya, size); + CHECK(initres,initres); + gsl_interp_accel *acc = gsl_interp_accel_alloc(); + if (NULL == acc) { gsl_spline_free(spline); ERROR(MEM); }; + + const int res = gsl_spline_eval_deriv2_e(spline, x, acc, y); + CHECK(res,res); + gsl_interp_accel_free(acc); + gsl_spline_free(spline); + OK +} + +int spline_eval_integ(const double xa[], const double ya[], unsigned int size, + double a, double b, int method, double *y) { + DEBUGMSG("spline_eval_integ"); + const gsl_interp_type *T; + switch (method) { + case 0: { T = gsl_interp_linear; break; } + case 1: { T = gsl_interp_polynomial; break; } + case 2: { T = gsl_interp_cspline; break; } + case 3: { T = gsl_interp_cspline_periodic; break; } + case 4: { T = gsl_interp_akima; break; } + case 5: { T = gsl_interp_akima_periodic; break; } + default: ERROR(BAD_CODE); + } + + gsl_spline *spline = gsl_spline_alloc(T, size); + if (NULL == spline) ERROR(MEM); + const int initres = gsl_spline_init(spline, xa, ya, size); + CHECK(initres,initres); + gsl_interp_accel *acc = gsl_interp_accel_alloc(); + if (NULL == acc) { gsl_spline_free(spline); ERROR(MEM); }; + + const int res = gsl_spline_eval_integ_e(spline, a, b, acc, y); + CHECK(res,res); + gsl_interp_accel_free(acc); + gsl_spline_free(spline); + OK +} + int integrate_qng(double f(double, void*), double a, double b, double aprec, double prec, double *result, double*error) { DEBUGMSG("integrate_qng"); @@ -363,7 +477,30 @@ OK } - +int siman(int seed, + gsl_siman_params_t *params, void *xp0, + double energy(void *), double metric(void *, void *), + void step(const gsl_rng *, void *, double), + void copy(void *, void *), void *copycons(void *), + void destroy(void *), void print(void *)) { + DEBUGMSG("siman"); + gsl_rng *gen = gsl_rng_alloc (gsl_rng_mt19937); + gsl_rng_set(gen, seed); + + // The simulated annealing routine doesn't indicate with a return + // code how things went -- there's little notion of convergence for + // a randomized minimizer on a potentially non-convex problem, and I + // suppose it doesn't detect egregious failures like malloc errors + // in the copy-constructor. + gsl_siman_solve(gen, xp0, + energy, step, + metric, print, + copy, copycons, + destroy, 0, *params); + + gsl_rng_free(gen); + OK +} // this version returns info about intermediate steps int minimize(int method, double f(int, double*), double tolsize, int maxit, diff --git a/src/Numeric/GSL/gsl-ode.c b/src/Numeric/GSL/gsl-ode.c --- a/src/Numeric/GSL/gsl-ode.c +++ b/src/Numeric/GSL/gsl-ode.c @@ -23,10 +23,11 @@ } -int ode(int method, double h, double eps_abs, double eps_rel, +int ode(int method, int control, double h, + double eps_abs, double eps_rel, double a_y, double a_dydt, int f(double, int, const double*, int, double*), int jac(double, int, const double*, int, int, double*), - KRVEC(xi), KRVEC(ts), RMAT(sol)) { + KRVEC(sc), KRVEC(xi), KRVEC(ts), RMAT(sol)) { const gsl_odeiv_step_type * T; @@ -46,9 +47,17 @@ } gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin); - gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel); gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin); + gsl_odeiv_control * c; + switch(control) { + case 0: { c = gsl_odeiv_control_standard_new + (eps_abs, eps_rel, a_y, a_dydt); break; } + case 1: { c = gsl_odeiv_control_scaled_new + (eps_abs, eps_rel, a_y, a_dydt, scp, scn); break; } + default: ERROR(BAD_CODE); + } + Tode P; P.f = f; P.j = jac; @@ -112,10 +121,11 @@ } -int ode(int method, double h, double eps_abs, double eps_rel, +int ode(int method, int control, double h, + double eps_abs, double eps_rel, double a_y, double a_dydt, int f(double, int, const double*, int, double*), int jac(double, int, const double*, int, int, double*), - KRVEC(xi), KRVEC(ts), RMAT(sol)) { + KRVEC(sc), KRVEC(xi), KRVEC(ts), RMAT(sol)) { const gsl_odeiv2_step_type * T; @@ -141,8 +151,15 @@ gsl_odeiv2_system sys = {odefunc, odejac, xin, &P}; - gsl_odeiv2_driver * d = - gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel); + gsl_odeiv2_driver * d; + + switch(control) { + case 0: { d = gsl_odeiv2_driver_alloc_standard_new + (&sys, T, h, eps_abs, eps_rel, a_y, a_dydt); break; } + case 1: { d = gsl_odeiv2_driver_alloc_scaled_new + (&sys, T, h, eps_abs, eps_rel, a_y, a_dydt, scp); break; } + default: ERROR(BAD_CODE); + } double t = tsp[0];