added odeSolve

11 files changed, 257 insertions, 8 deletions

diff --git a/CHANGES b/CHANGES
index 8dd8066..32d34c0 100644
--- a/CHANGES
+++ b/CHANGES
@@ -1,9 +1,11 @@
 0.8.3.0
 =======
 
-- Matrix arithmetic automatically replicates single row/column matrices.
+- odeSolve
 
-- latexFormat, dispcf for complex matrices
+- Matrix arithmetic automatically replicates matrix with single row/column
+
+- latexFormat, dispcf
 
 0.8.2.0
 =======
diff --git a/examples/HMatrixReal.hs b/examples/Real.hs
index 0edff10..b32a961 100644
--- a/examples/HMatrixReal.hs
+++ b/examples/Real.hs
@@ -1,7 +1,7 @@
 
 -- Alternative interface and utilities for creation of real arrays, useful to work in interactive mode.
 
-module HMatrixReal(
+module Real(
     module Numeric.LinearAlgebra,
     (<>), (*>), (<*), (<\>), (\>),
     vector,
diff --git a/examples/ode.hs b/examples/ode.hs
new file mode 100644
index 0000000..082c46c
--- /dev/null
+++ b/examples/ode.hs
@@ -0,0 +1,34 @@
+import Numeric.GSL.ODE
+import Numeric.LinearAlgebra
+import Graphics.Plot
+
+vanderpol mu = do
+    let xdot mu t [x,v] = [v, -x + mu * v * (1-x^2)]
+        ts = linspace 1000 (0,50)
+        sol = toColumns $ odeSolve (xdot mu) [1,0] ts
+    mplot (ts : sol)
+    mplot sol
+
+
+harmonic w d = do
+    let xdot w d t [x,v] = [v, a*x + b*v] where a = -w^2; b = -2*d*w
+        ts = linspace 100 (0,20)
+        sol = odeSolve (xdot w d) [1,0] ts
+    mplot (ts : toColumns sol)
+
+
+kepler v a = mplot (take 2 $ toColumns sol) where
+    xdot t [x,y,vx,vy] = [vx,vy,x*k,y*k]
+        where g=1
+              k=(-g)*(x*x+y*y)**(-1.5)
+    ts = linspace 100 (0,30)
+    sol = odeSolve xdot [4, 0, v * cos (a*degree), v * sin (a*degree)] ts
+    degree = pi/180
+
+
+main = do
+    vanderpol 2
+    harmonic 1 0
+    harmonic 1 0.1
+    kepler 0.3 60
+    kepler 0.4 70
diff --git a/hmatrix.cabal b/hmatrix.cabal
index 6f6c832..1080dcd 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -25,6 +25,7 @@ extra-source-files: examples/tests.hs
                     examples/integrate.hs
                     examples/minimize.hs
                     examples/root.hs
+                    examples/ode.hs
                     examples/pca1.hs
                     examples/pca2.hs
                     examples/pinv.hs
@@ -37,7 +38,7 @@ extra-source-files: examples/tests.hs
                     examples/error.hs
                     examples/devel/wrappers.hs
                     examples/devel/functions.c
-                    examples/HMatrixReal.hs
+                    examples/Real.hs
 
 extra-source-files: lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h,
                     lib/Numeric/LinearAlgebra/LAPACK/clapack.h
@@ -85,6 +86,7 @@ library
                         Numeric.GSL.Polynomials,
                         Numeric.GSL.Minimization,
                         Numeric.GSL.Root,
+                        Numeric.GSL.ODE,
                         Numeric.GSL.Vector,
                         Numeric.GSL.Special,
                         Numeric.GSL.Special.Gamma,
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index dced49f..e204a83 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -28,7 +28,7 @@ module Data.Packed.Matrix (
     extractRows,
     ident, diag, diagRect, takeDiag,
     liftMatrix, liftMatrix2, liftMatrix2Auto,
-    dispf, disps, dispcf, showComplex, latexFormat, format,
+    dispf, disps, dispcf, latexFormat, format,
     loadMatrix, saveMatrix, fromFile, fileDimensions,
     readMatrix, fromArray2D
 ) where
@@ -339,7 +339,7 @@ lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx
 isZero x = show x `elem` ["0.0","-0.0"]
 isOne  x = show x `elem` ["1.0","-1.0"]
 
--- | Pretty print a complex matrix with with at most n decimal digits.
+-- | Pretty print a complex matrix with at most n decimal digits.
 dispcf :: Int -> Matrix (Complex Double) -> String
 dispcf d m = sdims m ++ "\n" ++ format "  " (showComplex d) m
 
diff --git a/lib/Numeric/GSL.hs b/lib/Numeric/GSL.hs
index cc5d156..68430bc 100644
--- a/lib/Numeric/GSL.hs
+++ b/lib/Numeric/GSL.hs
@@ -19,6 +19,7 @@ module Numeric.GSL (
 , module Numeric.GSL.Polynomials
 , module Numeric.GSL.Minimization
 , module Numeric.GSL.Root
+, module Numeric.GSL.ODE
 , module Numeric.GSL.Special
 , module Complex
 , setErrorHandlerOff
@@ -31,6 +32,7 @@ import Numeric.GSL.Fourier
 import Numeric.GSL.Polynomials
 import Numeric.GSL.Minimization
 import Numeric.GSL.Root
+import Numeric.GSL.ODE
 import Complex
 
 
diff --git a/lib/Numeric/GSL/Internal.hs b/lib/Numeric/GSL/Internal.hs
index 834dfc2..37bcc1b 100644
--- a/lib/Numeric/GSL/Internal.hs
+++ b/lib/Numeric/GSL/Internal.hs
@@ -30,6 +30,9 @@ foreign import ccall "wrapper"
 foreign import ccall "wrapper"
     mkVecVecfun :: TVV -> IO (FunPtr TVV)
 
+foreign import ccall "wrapper"
+    mkDoubleVecVecfun :: (Double -> TVV) -> IO (FunPtr (Double -> TVV))
+
 aux_vTov :: (Vector Double -> Vector Double) -> TVV
 aux_vTov f n p nr r = g where
     V {fptr = pr} = f x
@@ -43,6 +46,9 @@ aux_vTov f n p nr r = g where
 foreign import ccall "wrapper"
     mkVecMatfun :: TVM -> IO (FunPtr TVM)
 
+foreign import ccall "wrapper"
+    mkDoubleVecMatfun :: (Double -> TVM) -> IO (FunPtr (Double -> TVM))
+
 aux_vTom :: (Vector Double -> Matrix Double) -> TVM
 aux_vTom f n p rr cr r = g where
     V {fptr = pr} = flatten $ f x
diff --git a/lib/Numeric/GSL/ODE.hs b/lib/Numeric/GSL/ODE.hs
new file mode 100644
index 0000000..3450b91
--- /dev/null
+++ b/lib/Numeric/GSL/ODE.hs
@@ -0,0 +1,107 @@
+{- |
+Module      :  Numeric.GSL.ODE
+Copyright   :  (c) Alberto Ruiz 2010
+License     :  GPL
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+Solution of ordinary differential equation (ODE) initial value problems.
+
+<http://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html>
+
+A simple example:
+
+@import Numeric.GSL
+import Numeric.LinearAlgebra
+import Graphics.Plot
+
+xdot t [x,v] = [v, -0.95*x - 0.1*v]
+
+ts = linspace 100 (0,20)
+
+sol = odeSolve xdot [10,0] ts
+
+main = mplot (ts : toColumns sol)@
+
+-}
+-----------------------------------------------------------------------------
+
+module Numeric.GSL.ODE (
+    odeSolve, odeSolveV, ODEMethod(..)
+) where
+
+import Data.Packed.Internal
+import Data.Packed.Matrix
+import Foreign
+import Foreign.C.Types(CInt)
+import Numeric.GSL.Internal
+
+-------------------------------------------------------------------------
+
+-- | Stepping functions
+data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method.
+               | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use 'RKf45'.
+               | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator.
+               | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method.
+               | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method.
+               | RK2imp -- ^ Implicit 2nd order Runge-Kutta at Gaussian points.
+               | RK4imp -- ^ Implicit 4th order Runge-Kutta at Gaussian points.
+               | BSimp -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian.
+               | Gear1 -- ^ M=1 implicit Gear method.
+               | Gear2 -- ^ M=2 implicit Gear method.
+               deriving (Enum,Eq,Show,Bounded)
+
+-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists.
+odeSolve
+    :: (Double -> [Double] -> [Double])        -- ^ xdot(t,x)
+    -> [Double]        -- ^ initial conditions
+    -> Vector Double   -- ^ desired solution times
+    -> Matrix Double   -- ^ solution
+odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) Nothing (fromList xi) ts
+    where hi = (ts@>1 - ts@>0)/100
+          epsAbs = 1.49012e-08
+          epsRel = 1.49012e-08
+          l2v f = \t -> fromList  . f t . toList
+          l2m f = \t -> fromLists . f t . toList
+
+-- | Evolution of the system with adaptive step-size 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)
+    -> Maybe (Double -> Vector Double -> Matrix Double)   -- ^ optional jacobian
+    -> Vector Double     -- ^ initial conditions
+    -> Vector Double     -- ^ desired solution times
+    -> Matrix Double     -- ^ solution
+odeSolveV method h epsAbs epsRel f mbjac 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 <- withVector xiv $ \xiv' ->
+            withVector ts $ \ts' ->
+             createMIO (dim ts) n
+              (ode_c (fi (fromEnum method)) h epsAbs epsRel fp jp // xiv' // ts' )
+              "ode"
+    freeHaskellFunPtr fp
+    return sol
+
+foreign import ccall "ode"
+    ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM
+
+-------------------------------------------------------
+
+checkdim1 n v
+    | dim v == n = v
+    | otherwise = error $ "Error: "++ show n
+                        ++ " components expected in the result of the function supplied to odeSolve"
+
+checkdim2 n m
+    | rows m == n && cols m == n = m
+    | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
+                        ++ " Jacobian expected in odeSolve"
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 0c71ca1..0fcde1c 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -22,6 +22,7 @@
 #include <gsl/gsl_complex_math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_odeiv.h>
 #include <string.h>
 #include <stdio.h>
 
@@ -507,7 +508,7 @@ double f_aux_min(const gsl_vector*x, void *pars) {
 
 
 void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) {
-    Tfdf * fdf = ((Tfdf*) pars);  
+    Tfdf * fdf = ((Tfdf*) pars);
     double* p = (double*)calloc(x->size,sizeof(double));
     double* q = (double*)calloc(g->size,sizeof(double));
     int k;
@@ -797,3 +798,89 @@ int random_vector(int seed, int code, RVEC(r)) {
     }
 }
 #undef RAN
+
+//////////////////////////////////////////////////////
+//                        ODE
+
+typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;
+
+int odefunc (double t, const double y[], double f[], void *params) { 
+    Tode * P = (Tode*) params;
+    (P->f)(t,P->n,y,P->n,f);
+    return GSL_SUCCESS;
+}
+
+int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
+     Tode * P = ((Tode*) params);
+     (P->j)(t,P->n,y,P->n,P->n,dfdy);
+     int j;
+     for (j=0; j< P->n; j++)
+        dfdt[j] = 0.0;
+     return GSL_SUCCESS;
+}
+
+
+int ode(int method, double h, double eps_abs, double eps_rel,
+        int f(double, int, const double*, int, double*),
+        int jac(double, int, const double*, int, int, double*),
+        KRVEC(xi), KRVEC(ts), RMAT(sol)) {
+
+    const gsl_odeiv_step_type * T;
+
+    switch(method) {
+        case 0 : {T = gsl_odeiv_step_rk2; break; }
+        case 1 : {T = gsl_odeiv_step_rk4; break; }
+        case 2 : {T = gsl_odeiv_step_rkf45; break; }
+        case 3 : {T = gsl_odeiv_step_rkck; break; }
+        case 4 : {T = gsl_odeiv_step_rk8pd; break; }
+        case 5 : {T = gsl_odeiv_step_rk2imp; break; }
+        case 6 : {T = gsl_odeiv_step_rk4imp; break; }
+        case 7 : {T = gsl_odeiv_step_bsimp; break; }
+        case 8 : {T = gsl_odeiv_step_gear1; break; }
+        case 9 : {T = gsl_odeiv_step_gear2; break; }
+        default: ERROR(BAD_CODE);
+    }
+
+
+    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);
+
+    Tode P;
+    P.f = f;
+    P.j = jac;
+    P.n = xin;
+
+    gsl_odeiv_system sys = {odefunc, odejac, xin, &P};
+
+    double t = tsp[0];
+
+    double* y = (double*)calloc(xin,sizeof(double));
+    int i,j;
+    for(i=0; i< xin; i++) {
+        y[i] = xip[i];
+        solp[i] = xip[i];
+    }
+
+       for (i = 1; i < tsn ; i++)
+         {
+           double ti = tsp[i];
+           while (t < ti)
+             {
+               gsl_odeiv_evolve_apply (e, c, s,
+                                       &sys,
+                                       &t, ti, &h,
+                                       y);
+               // if (h < hmin) h = hmin;
+             }
+           for(j=0; j<xin; j++) {
+               solp[i*xin + j] = y[j];
+           }
+         }
+
+    free(y);
+    gsl_odeiv_evolve_free (e);
+    gsl_odeiv_control_free (c);
+    gsl_odeiv_step_free (s);
+    return 0;
+}
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index e6b26a0..46c1804 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -135,6 +135,14 @@ rootFindingTest = TestList [ utest "root Hybrids" (fst sol1 ~~ [1,1])
 
 ---------------------------------------------------------------------
 
+odeTest = utest "ode" (last (toLists sol) ~~ [-1.7588880332411019, 8.364348908711941e-2])
+    where sol = odeSolveV RK8pd 1E-6 1E-6 0 (l2v $ vanderpol 10) Nothing (fromList [1,0]) ts
+          ts = linspace 101 (0,100)
+          l2v f = \t -> fromList  . f t . toList
+          vanderpol mu _t [x,y] = [y, -x + mu * y * (1-x^2) ]
+
+---------------------------------------------------------------------
+
 randomTestGaussian = c :~1~: snd (meanCov dat) where
     a = (3><3) [1,2,3,
                 2,4,0,
@@ -280,6 +288,7 @@ runTests n = do
         , utest "rank" $  rank ((2><3)[1,0,0,1,6*eps,0]) == 1
                        && rank ((2><3)[1,0,0,1,7*eps,0]) == 2
         , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM)
+        , odeTest
         ]
     return ()
 
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
index 405ce64..42c6ce2 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -148,7 +148,7 @@ svdProp4 m' = m |~| u <> real (diag s) <> trans v
            && orthonormal u && orthonormal v
            && (dim s == r || r == 0 && dim s == 1)
     where (u,s,v) = compactSVD m
-          m = m' <-> m'
+          m = fromBlocks [[m'],[m']]
           r = rank m'
 
 svdProp5a m = and (map (s1|~|) [s2,s3,s4,s5,s6]) where