hmatrix-gsl 0.16.0.3 → 0.17.0.0

raw patch · 18 files changed

+874/−126 lines, 18 filesdep ~hmatrixPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: hmatrix

API changes (from Hackage documentation)

Files

• hmatrix-gsl.cabal +10/−3
• src/Graphics/Plot.hs +2/−2
• src/Numeric/GSL.hs +2/−0
• src/Numeric/GSL/Fitting.hs +10/−8
• src/Numeric/GSL/Fourier.hs +6/−4
• src/Numeric/GSL/IO.hs +1/−1
• src/Numeric/GSL/Internal.hs +19/−8
• src/Numeric/GSL/Interpolation.hs +284/−0
• src/Numeric/GSL/LinearAlgebra.hs +7/−7
• src/Numeric/GSL/Minimization.hs +11/−8
• src/Numeric/GSL/ODE.hs +85/−49
• src/Numeric/GSL/Polynomials.hs +4/−5
• src/Numeric/GSL/Random.hs +9/−7
• src/Numeric/GSL/Root.hs +10/−8
• src/Numeric/GSL/SimulatedAnnealing.hs +245/−0
• src/Numeric/GSL/Vector.hs +7/−8
• src/Numeric/GSL/gsl-aux.c +138/−1
• src/Numeric/GSL/gsl-ode.c +24/−7

@@ -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

@@ -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.

@@ -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  

@@ -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)  

@@ -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+

@@ -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) 

@@ -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 (#) #-} 

@@ -0,0 +1,284 @@+{- |+Module      :  Numeric.GSL.Interpolation+Copyright   :  (c) Matthew Peddie 2015+License     :  GPL+Maintainer  :  Alberto Ruiz+Stability   :  provisional++Interpolation routines.++<https://www.gnu.org/software/gsl/manual/html_node/Interpolation.html#Interpolation>++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

@@ -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

@@ -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"

@@ -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

@@ -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" 

@@ -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

@@ -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" 

@@ -0,0 +1,245 @@+{- |+Module      :  Numeric.GSL.Interpolation+Copyright   :  (c) Matthew Peddie 2015+License     :  GPL+Maintainer  :  Alberto Ruiz+Stability   :  provisional++Simulated annealing routines.++<https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing.html#Simulated-Annealing>++Here is a translation of the simple example given in+<https://www.gnu.org/software/gsl/manual/html_node/Trivial-example.html#Trivial-example the GSL manual>:++> 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+-- <https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing-functions.html#Simulated-Annealing-functions the GSL manual>,+-- which controls the simulated annealing algorithm.+--+-- The annealing process is parameterized by the Boltzmann+-- distribution and the /cooling schedule/.  For more details, see+-- <https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing-algorithm.html#Simulated-Annealing-algorithm the relevant section of the manual>.+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

@@ -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

@@ -34,7 +34,10 @@ #include <gsl/gsl_rng.h> #include <gsl/gsl_randist.h> #include <gsl/gsl_roots.h>+#include <gsl/gsl_spline.h> #include <gsl/gsl_multifit_nlin.h>+#include <gsl/gsl_siman.h>+ #include <string.h> #include <stdio.h> @@ -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, 

@@ -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];