Merge pull request #2 from owlbarn/first-class-ode

First class ode

.travis.yml

@@ -6,6 +6,7 @@ install: wget https://raw.githubusercontent.com/ocaml/ocaml-travisci-skeleton/ma
 script: bash -ex ./.travis-docker.sh
 env:
   global:
-  - PINS="ode:."
+  - PINS="owl-ode:. owl-ode-sundials:."
   matrix:
-  - PACKAGE="ode" DISTRO="debian-stable" OCAML_VERSION="4.07"
+  - PACKAGE="owl-ode" DISTRO="debian-stable" OCAML_VERSION="4.07"
+  - PACKAGE="owl-ode-sundials" DISTRO="debian-stable" OCAML_VERSION="4.07"

examples/damped.ml

@@ -1,4 +1,7 @@
-let damped_noforcing a xs ps _ : Owl.Mat.mat =
+open Owl_ode
+open Owl_ode.Types
+
+let damped_noforcing a (xs, ps) _ : Owl.Mat.mat =
   Owl.Mat.(
     xs *$ (-1.0) + ps *$ (-1.0*.a)
   )
@@ -25,12 +28,9 @@ let () =
   let p0 = Owl.Mat.of_array [|0.75|] 1 1 in
   let t0, duration = 0.0, 15.0 in
   let f = damped_noforcing a in
-  let open Ode.SymplecticSolver in
-  let tspec = Ode.T1 {t0; duration; dt} in
-  let problem = {f; x0; p0} in
-  let ndsolve algo = odeint ~algo ~problem ~tspec in
-  let t, sol1, _ = ndsolve leapfrog () in
-  let _, sol2, _ = ndsolve ruth3 () in
-  let _, sol3, _ = ndsolve ruth4 () in
+  let tspec = T1 {t0; duration; dt} in
+  let t, sol1, _ = Ode.odeint (module Symplectic.Leapfrog) f (x0,p0) tspec () in
+  let _, sol2, _ = Ode.odeint (module Symplectic.Ruth3) f (x0,p0) tspec () in
+  let _, sol3, _ = Ode.odeint (module Symplectic.Symplectic_Euler) f (x0,p0) tspec () in
   (* XXX: I'd prefer t to be already an Owl array as well *)
   plot_sol "damped.png" (Owl.Mat.of_array t 1 (Array.length t)) sol1 sol2 sol3;

examples/dune

@@ -1,11 +1,11 @@
 (executable
  (name van_der_pol)
- (libraries ode)
+ (libraries owl owl-ode-sundials)
  (modules van_der_pol)
 )
 
 (executable
  (name damped)
- (libraries ode)
+ (libraries owl owl-ode)
  (modules damped)
 )

examples/van_der_pol.ml

@@ -1,4 +1,6 @@
 open Owl
+open Owl_ode
+open Owl_ode.Types
 
 let () = Printexc.record_backtrace true
 
@@ -14,12 +16,26 @@ let f =
 
 let y0 = Mat.of_array [|0.02; 0.03|] 2 1
 
+(* use Ode provided cvode integrator *)
 let () = 
-  let open Ode.OdeSolver in
-  let problem = {f; y0} in
-  let tspec = Ode.T1 {t0=0.0; dt=1E-2; duration=30.0} in
-  let ts, ys = Ode.OdeSolver.odeint ~algo:(cvode ()) ~problem ~tspec () in
+  let tspec = T1 {t0=0.0; dt=1E-2; duration=30.0} in
+  let ts, ys = Ode.odeint (module Owl_ode_sundials.Owl_Cvode) f y0 tspec () in
   (* save ts and ys *)
   let ts = [| ts |] |> Mat.of_arrays |> Mat.transpose in
   let ys = ys |> Mat.transpose in
   Mat.save_txt Mat.(ts @|| ys) "van_der_pol_dynamics.txt"
+
+(* create our own cvode integrator *)
+module Custom_Cvode = struct
+  type t = Mat.mat
+  type output = float array * Mat.mat
+  let solve = Owl_ode_sundials.cvode ~stiff:false ~relative_tol:1E-3 ()
+end
+
+let () = 
+  let tspec = T1 {t0=0.0; dt=1E-2; duration=30.0} in
+  let ts, ys = Ode.odeint (module Custom_Cvode) f y0 tspec () in
+  (* save ts and ys *)
+  let ts = [| ts |] |> Mat.of_arrays |> Mat.transpose in
+  let ys = ys |> Mat.transpose in
+  Mat.save_txt Mat.(ts @|| ys) "van_der_pol_dynamics_custom.txt"

owl-ode-sundials.opam

@@ -0,0 +1,20 @@
+opam-version: "2.0"
+maintainer: "owlbarn"
+authors: [ "Marcello Seri" "Ta-Chu Calvin Kao" ]
+license: "MIT"
+homepage: "https://github.com/owlbarn/owl_sundials"
+dev-repo: "git+https://github.com/owlbarn/owl_sundials.git"
+bug-reports: "https://github.com/owlbarn/owl_sundials/issues"
+doc: "https://owlbarn.github.io/owl_sundials/ode"
+build: [
+  ["dune" "build" "-p" name "-j" jobs]
+  ["dune" "exec" "examples/van_der_pol.exe"] {with-test}
+]
+depends: [
+  "ocaml" {>= "4.02"}
+  "owl" {>= "0.4.0"}
+  "dune" {build & >= "1.1.0"}
+  "owl-ode" (*{>="0.0.4"}*)
+  "sundialsml"
+]
+synopsis: "Owl's ODE solvers, interface with SundialsML"

ode.opam → owl-ode.opam

@@ -8,13 +8,11 @@ bug-reports: "https://github.com/owlbarn/owl_sundials/issues"
 doc: "https://owlbarn.github.io/owl_sundials/ode"
 build: [
   ["dune" "build" "-p" name "-j" jobs]
-  ["dune" "exec" "examples/van_der_pol.exe"] {with-test}
   ["dune" "exec" "examples/damped.exe"] {with-test}
 ]
 depends: [
   "ocaml" {>= "4.02"}
   "owl" {>= "0.4.0"}
   "dune" {build & >= "1.1.0"}
-  "sundialsml"
 ]
 synopsis: "Owl's ODE solvers"

src/common.ml → src/ode/common.ml

@@ -2,9 +2,6 @@ open Owl
 (* TODO: update implementations of multiple order RK on the line of
  * symplectic.ml *)
 
-type timespan = float * float
-(** Representation of a time span. *)  
-
 let steps t0 t1 dt =
   (* NOTE: switched Float.floor to Maths.floor; 
    * Float module seems not to be only supported in ocaml 4.07.0 *)
@@ -44,8 +41,6 @@ let integrate ~step ~tspan:(t0, t1) ~dt y0 =
   !ts |> List.rev |> Array.of_list,
   ys
 
-
-
 let symplectic_integrate ~step ~tspan:(t0, t1) ~dt x0 p0 =
   assert ((Mat.shape x0)=(Mat.shape p0));
   let major, n = get_major x0 in
@@ -73,4 +68,3 @@ let symplectic_integrate ~step ~tspan:(t0, t1) ~dt x0 p0 =
   !ts |> List.rev |> Array.of_list,
   xs, ps 
 
-

src/ode/dune

@@ -0,0 +1,5 @@
+(library
+ (name owl_ode)
+ (public_name owl-ode)
+ (libraries owl)
+)

src/ode/native.ml

@@ -0,0 +1,41 @@
+open Owl
+open Types
+
+type f_t = Mat.mat -> float -> Mat.mat
+
+let euler_s ~(f:f_t) ~dt = fun y0 t0 ->
+  let y = Mat.(y0 + (f y0 t0) *$ dt) in
+  let t = t0 +. dt in
+  y, t
+
+let rk4_s ~(f:f_t) ~dt = fun y0 t0 ->
+  let k1 = Mat.(dt $* (f y0 t0)) in
+  let k2 = Mat.(dt $* (f (y0 + k1 *$ 0.5) (t0 +. 0.5 *. dt))) in
+  let k3 = Mat.(dt $* (f (y0 + k2 *$ 0.5) (t0 +. 0.5 *. dt))) in
+  let k4 = Mat.(dt $* (f (y0 + k3) (t0 +. dt))) in
+  let dy = Mat.((k1 + k2 + k3 + k4) /$ 6.) in
+  let y = Mat.(y0 + dy) in
+  let t = t0 +. dt in
+  y, t
+
+let prepare step f y0 tspec () =
+  let tspan, dt = match tspec with
+    | T1 {t0; duration; dt} -> (t0, t0+.duration), dt
+    | T2 {tspan; dt} -> tspan, dt 
+    | T3 _ -> raise Owl_exception.NOT_IMPLEMENTED 
+  in
+  let step = step ~f ~dt in
+  Common.integrate ~step ~tspan ~dt y0
+
+
+module Euler = struct
+  type t = Mat.mat
+  type output = float array * Mat.mat
+  let solve = prepare euler_s
+end
+
+module RK4 = struct
+  type t = Mat.mat
+  type output = float array * Mat.mat
+  let solve = prepare rk4_s
+end