Using hsc, semi-auto-generating documentation

7 files changed, 303 insertions, 86 deletions

diff --git a/packages/sundials/hmatrix-sundials.cabal b/packages/sundials/hmatrix-sundials.cabal
index 762537e..dbe50e0 100644
--- a/packages/sundials/hmatrix-sundials.cabal
+++ b/packages/sundials/hmatrix-sundials.cabal
@@ -16,8 +16,6 @@ build-type:          Simple
 extra-source-files:  ChangeLog.md
 cabal-version:       >=1.10
 
-extra-source-files:  src/helpers.c, src/helpers.h
-
 
 library
   build-depends:       base >=4.10 && <4.11,
@@ -26,15 +24,19 @@ library
                        template-haskell >=2.12 && <2.13,
                        containers >=0.5 && <0.6,
                        hmatrix>=0.18
-  other-extensions:    QuasiQuotes, TemplateHaskell, MultiWayIf, OverloadedStrings
+  extra-libraries:     sundials_arkode
+  other-extensions:    QuasiQuotes
   hs-source-dirs:      src
   exposed-modules:     Numeric.Sundials.Arkode.ODE
-  other-modules:       Types
+  other-modules:       Types,
+                       Bar
 
 executable sundials
   main-is:             Main.hs
-  other-modules:       Types, Numeric.Sundials.Arkode.ODE
-  other-extensions:    QuasiQuotes, TemplateHaskell, MultiWayIf, OverloadedStrings
+  other-modules:       Types,
+                       Numeric.Sundials.Arkode.ODE,
+                       Bar
+  other-extensions:    QuasiQuotes
   build-depends:       base >=4.10 && <4.11,
                        inline-c >=0.6 && <0.7,
                        vector >=0.12 && <0.13,
@@ -43,7 +45,9 @@ executable sundials
                        hmatrix>=0.18,
                        plots,
                        diagrams-lib,
-                       diagrams-rasterific
+                       diagrams-rasterific,
+                       lens,
+                       pretty
   hs-source-dirs:      src
   default-language:    Haskell2010
   extra-libraries:     sundials_arkode
diff --git a/packages/sundials/src/Bar.hsc b/packages/sundials/src/Bar.hsc
index 434c4d4..7db0d4a 100644
--- a/packages/sundials/src/Bar.hsc
+++ b/packages/sundials/src/Bar.hsc
@@ -1,6 +1,6 @@
 {-# LANGUAGE RecordWildCards #-}
 
-module Example where
+module Bar where
 
 import Foreign
 import Foreign.C.Types
@@ -8,6 +8,7 @@ import Foreign.C.String
 
 #include "/Users/dom/sundials/include/sundials/sundials_nvector.h"
 #include "/Users/dom/sundials/include/nvector/nvector_serial.h"
+#include "/Users/dom/sundials/include/arkode/arkode.h"
 
 #def typedef struct _generic_N_Vector BarType;
 #def typedef struct _N_VectorContent_Serial BazType;
@@ -19,6 +20,9 @@ getContentPtr ptr = (#peek BarType, content) ptr
 getData :: Storable a => Ptr b -> IO a
 getData ptr = (#peek BazType, data) ptr
 
+arkSMax :: Int
+arkSMax = #const ARK_S_MAX
+
 
   
 
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index 5e51372..895e610 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -11,6 +11,9 @@ import           Diagrams.Backend.Rasterific
 import           Control.Lens
 import           Data.List (zip4)
 
+import           Text.PrettyPrint.HughesPJClass
+import           Data.List (intercalate)
+
 
 brusselator _t x = [ a - (w + 1) * u + v * u^2
                    , w * u - v * u^2
@@ -43,8 +46,37 @@ kSaxis :: [(Double, Double)] -> P.Axis B D.V2 Double
 kSaxis xs = P.r2Axis &~ do
   P.linePlot' xs
 
+butcherTableauTex :: (Show a, Element a) => Matrix a -> String
+butcherTableauTex m = render $
+                    vcat [ text ("\n\\begin{array}{c|" ++ (concat $ replicate n "c") ++ "}")
+                         , us
+                         , text "\\end{array}"
+                         ]
+  where
+    n = rows m
+    rs = toLists m
+    ss = map (\r -> intercalate " & " $ map show r) rs
+    ts = zipWith (\n r -> "c_" ++ show n ++ " & " ++ r) [1..n] ss
+    us = vcat $ map (\r -> text r <+> text "\\\\") ts
+
 main :: IO ()
 main = do
+  -- $$
+  -- \begin{array}{c|cccc}
+  -- c_1    & a_{11} & a_{12}& \dots & a_{1s}\\
+  -- c_2    & a_{21} & a_{22}& \dots & a_{2s}\\
+  -- \vdots & \vdots & \vdots& \ddots& \vdots\\
+  -- c_s    & a_{s1} & a_{s2}& \dots & a_{ss} \\
+  -- \hline
+  --        & b_1    & b_2   & \dots & b_s\\
+  --        & b^*_1  & b^*_2 & \dots & b^*_s\\
+  -- \end{array}
+  -- $$
+
+  let res = btGet
+  putStrLn $ show res
+  putStrLn $ butcherTableauTex res
+  
   let res = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
   putStrLn $ show res
   renderRasterific "diagrams/brusselator.png"
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
index f432951..76ed61b 100644
--- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -6,8 +6,25 @@
 {-# LANGUAGE OverloadedStrings #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 
+-- |
+-- Module: Numeric.Sundials.ARKode
+--
+-- Blah
+--
+-- \[
+-- \begin{array}{c|cccc}
+-- c_1 & 0.0 & 0.0 & 0.0 & 0.0 \\
+-- c_2 & 0.4358665215 & 0.4358665215 & 0.0 & 0.0 \\
+-- c_3 & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\
+-- c_4 & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\
+-- \end{array}
+-- \]
+--
 module Numeric.Sundials.Arkode.ODE ( solveOde
                                    , odeSolve
+                                   , getButcherTable
+                                   , getBT
+                                   , btGet
                                    ) where
 
 import qualified Language.C.Inline as C
@@ -28,9 +45,10 @@ import           System.IO.Unsafe (unsafePerformIO)
 
 import           Numeric.LinearAlgebra.Devel (createVector)
 
-import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, (><))
+import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, (><), subMatrix)
 
 import qualified Types as T
+import qualified Bar as B
 
 C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx)
 
@@ -83,16 +101,16 @@ vectorToC vec len ptr = do
   V.copy (VM.unsafeFromForeignPtr0 ptr' len) vec
 
 data SundialsDiagnostics = SundialsDiagnostics {
-    aRKodeGetNumSteps :: Int
-  , aRKodeGetNumStepAttempts :: Int
-  , aRKodeGetNumRhsEvals_fe :: Int
-  , aRKodeGetNumRhsEvals_fi :: Int
-  , aRKodeGetNumLinSolvSetups :: Int
-  , aRKodeGetNumErrTestFails :: Int
-  , aRKodeGetNumNonlinSolvIters :: Int
+    aRKodeGetNumSteps               :: Int
+  , aRKodeGetNumStepAttempts        :: Int
+  , aRKodeGetNumRhsEvals_fe         :: Int
+  , aRKodeGetNumRhsEvals_fi         :: Int
+  , aRKodeGetNumLinSolvSetups       :: Int
+  , aRKodeGetNumErrTestFails        :: Int
+  , aRKodeGetNumNonlinSolvIters     :: Int
   , aRKodeGetNumNonlinSolvConvFails :: Int
-  , aRKDlsGetNumJacEvals :: Int
-  , aRKDlsGetNumRhsEvals :: Int
+  , aRKDlsGetNumJacEvals            :: Int
+  , aRKDlsGetNumRhsEvals            :: Int
   } deriving Show
 
 odeSolve :: (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
@@ -312,3 +330,107 @@ solveOdeC fun f0 ts = unsafePerformIO $ do
     else do
       return $ Left res
 
+btGet :: Matrix Double
+btGet = case getBT of
+          Left c -> error $ show c -- FIXME
+          Right (v, sqp) -> subMatrix (0, 0) (4, 4) $ (B.arkSMax >< B.arkSMax) (V.toList v)
+
+getBT :: Either Int (V.Vector Double, V.Vector Int)
+getBT = case getButcherTable of
+          Left c -> Left $ fromIntegral c
+          Right (v, sqp) -> Right $ (coerce v, V.map fromIntegral sqp)
+
+getButcherTable :: Either CInt ((V.Vector CDouble), V.Vector CInt)
+getButcherTable = unsafePerformIO $ do
+  -- arkode seems to want an ODE in order to set and then get the
+  -- Butcher tableau so here's one to keep it happy
+  let fun :: CDouble -> V.Vector CDouble -> V.Vector CDouble
+      fun t ys = V.fromList [ ys V.! 0 ]
+      f0       = V.fromList [ 1.0 ]
+      ts       = V.fromList [ 0.0 ]
+      dim = V.length f0
+      nEq :: CLong
+      nEq = fromIntegral dim
+
+  -- FIXME: I believe these gets taken from the ghc heap and so should
+  -- be subject to garbage collection.
+  btSQP :: V.Vector CInt <- createVector 3
+  btSQPMut <- V.thaw btSQP
+  btAs :: V.Vector CDouble <- createVector (B.arkSMax * B.arkSMax)
+  btAsMut <- V.thaw btAs
+  -- We need the types that sundials expects. These are tied together
+  -- in 'Types'. FIXME: The Haskell type is currently empty!
+  let funIO :: CDouble -> Ptr T.BarType -> Ptr T.BarType -> Ptr () -> IO CInt
+      funIO x y f _ptr = do
+        -- Convert the pointer we get from C (y) to a vector, and then
+        -- apply the user-supplied function.
+        fImm <- fun x <$> getDataFromContents dim y
+        -- Fill in the provided pointer with the resulting vector.
+        putDataInContents fImm dim f
+        -- I don't understand what this comment means
+        -- Unsafe since the function will be called many times.
+        [CU.exp| int{ 0 } |]
+  res <- [C.block| int {
+                         /* general problem variables */
+                         int flag;                       /* reusable error-checking flag */
+                         N_Vector y = NULL;              /* empty vector for storing solution */
+                         void *arkode_mem = NULL;        /* empty ARKode memory structure */
+
+                         /* general problem parameters */
+                         /* initial time */
+                         realtype T0 = RCONST(($vec-ptr:(double *ts))[0]);
+                         /* number of dependent vars. */
+                         sunindextype NEQ = $(sunindextype nEq);
+
+                         /* Initialize data structures */
+                         y = N_VNew_Serial(NEQ);         /* Create serial vector for solution */
+                         if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
+                         /* Specify initial condition */
+                         int i, j;
+                         for (i = 0; i < NEQ; i++) {
+                           NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i];
+                         };
+                         arkode_mem = ARKodeCreate(); /* Create the solver memory */
+                         if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
+
+                         flag = ARKodeInit(arkode_mem, NULL, $fun:(int (* funIO) (double t, BarType y[], BarType dydt[], void * params)), T0, y);
+                         if (check_flag(&flag, "ARKodeInit", 1)) return 1;
+
+                         flag = ARKodeSetIRKTableNum(arkode_mem, KVAERNO_4_2_3);
+                         if (check_flag(&flag, "ARKode", 1)) return 1;
+
+                         int s, q, p;
+                         realtype *ai = (realtype *)malloc(ARK_S_MAX * ARK_S_MAX * sizeof(realtype));
+                         realtype *ae = (realtype *)malloc(ARK_S_MAX * ARK_S_MAX * sizeof(realtype));
+                         realtype *ci = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
+                         realtype *ce = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
+                         realtype *bi = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
+                         realtype *be = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
+                         realtype *b2i = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
+                         realtype *b2e = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
+                         flag = ARKodeGetCurrentButcherTables(arkode_mem, &s, &q, &p, ai, ae, ci, ce, bi, be, b2i, b2e);
+                         if (check_flag(&flag, "ARKode", 1)) return 1;
+                         $vec-ptr:(int *btSQPMut)[0] = s;
+                         $vec-ptr:(int *btSQPMut)[1] = q;
+                         $vec-ptr:(int *btSQPMut)[2] = p;
+                         for (i = 0; i < s; i++) {
+                           for (j = 0; j < s; j++) {
+                             /* FIXME: double should be realtype */
+                             ($vec-ptr:(double *btAsMut))[i * ARK_S_MAX + j] = ai[i * ARK_S_MAX + j];
+                           }
+                         }
+
+                         /* Clean up and return */
+                         N_VDestroy(y);            /* Free y vector */
+                         ARKodeFree(&arkode_mem);  /* Free integrator memory */
+ 
+                         return flag;
+                       } |]
+  if res == 0
+    then do
+      x <- V.freeze btAsMut
+      y <- V.freeze btSQPMut
+      return $ Right (x, y)
+    else do
+      return $ Left res
+
diff --git a/packages/sundials/stack.yaml b/packages/sundials/stack.yaml
deleted file mode 100644
index 9c6b17c..0000000
--- a/packages/sundials/stack.yaml
+++ /dev/null
@@ -1,66 +0,0 @@
-# This file was automatically generated by 'stack init'
-#
-# Some commonly used options have been documented as comments in this file.
-# For advanced use and comprehensive documentation of the format, please see:
-# https://docs.haskellstack.org/en/stable/yaml_configuration/
-
-# Resolver to choose a 'specific' stackage snapshot or a compiler version.
-# A snapshot resolver dictates the compiler version and the set of packages
-# to be used for project dependencies. For example:
-#
-# resolver: lts-3.5
-# resolver: nightly-2015-09-21
-# resolver: ghc-7.10.2
-# resolver: ghcjs-0.1.0_ghc-7.10.2
-# resolver:
-#  name: custom-snapshot
-#  location: "./custom-snapshot.yaml"
-resolver: lts-10.9
-
-# User packages to be built.
-# Various formats can be used as shown in the example below.
-#
-# packages:
-# - some-directory
-# - https://example.com/foo/bar/baz-0.0.2.tar.gz
-# - location:
-#    git: https://github.com/commercialhaskell/stack.git
-#    commit: e7b331f14bcffb8367cd58fbfc8b40ec7642100a
-# - location: https://github.com/commercialhaskell/stack/commit/e7b331f14bcffb8367cd58fbfc8b40ec7642100a
-#   extra-dep: true
-#  subdirs:
-#  - auto-update
-#  - wai
-#
-# A package marked 'extra-dep: true' will only be built if demanded by a
-# non-dependency (i.e. a user package), and its test suites and benchmarks
-# will not be run. This is useful for tweaking upstream packages.
-packages:
-- .
-# Dependency packages to be pulled from upstream that are not in the resolver
-# (e.g., acme-missiles-0.3)
-# extra-deps: []
-
-# Override default flag values for local packages and extra-deps
-# flags: {}
-
-# Extra package databases containing global packages
-# extra-package-dbs: []
-
-# Control whether we use the GHC we find on the path
-# system-ghc: true
-#
-# Require a specific version of stack, using version ranges
-# require-stack-version: -any # Default
-# require-stack-version: ">=1.6"
-#
-# Override the architecture used by stack, especially useful on Windows
-# arch: i386
-# arch: x86_64
-#
-# Extra directories used by stack for building
-# extra-include-dirs: [/path/to/dir]
-# extra-lib-dirs: [/path/to/dir]
-#
-# Allow a newer minor version of GHC than the snapshot specifies
-# compiler-check: newer-minor
\ No newline at end of file
diff --git a/shell.nix b/shell.nix
new file mode 100644
index 0000000..9566524
--- /dev/null
+++ b/shell.nix
@@ -0,0 +1,115 @@
+# From Mathieu's PR
+
+# {ghc}:
+
+# with import <nixpkgs> {};
+
+# let
+#   gfortranlibdir = "${gfortran7.cc.lib}";
+#   # XXX Workaround https://ghc.haskell.org/trac/ghc/ticket/11042.
+#   libHack = if stdenv.isDarwin then {
+#       DYLD_LIBRARY_PATH = [gfortranlibdir];
+#     } else {
+#       LD_LIBRARY_PATH = [gfortranlibdir];
+#     };
+# in
+# haskell.lib.buildStackProject ({
+#   name = "haskell-sundials";
+#   buildInputs = [blas liblapack sundials zlib];
+#   inherit ghc;
+#   # XXX Workaround https://ghc.haskell.org/trac/ghc/ticket/11042.
+#   extraArgs = ["--extra-lib-dirs=${gfortranlibdir}"];
+# } // libHack)
+  
+
+{ nixpkgs ? import <nixpkgs>
+  {}
+, compiler ? "ghc822"
+, doBenchmark ? false }:
+
+let
+
+  inherit (nixpkgs) pkgs;
+
+f = { mkDerivation, ad, base, diagrams-lib, diagrams-rasterific, gcc
+    , hmatrix, hmatrix-gsl, inline-c, plots, pretty, stdenv, sundials, vector }:
+
+mkDerivation {
+  pname = "haskell-sundials";
+  version = "0.1.0.0";
+  src = ./.;
+  isLibrary = false;
+  isExecutable = true;
+  executableHaskellDepends = [
+    ad
+    base
+    diagrams-lib
+    diagrams-rasterific
+    gcc
+    hmatrix
+    hmatrix-gsl
+    inline-c
+    plots
+    pretty
+    vector
+  ] ++ (if pkgs.stdenv.isDarwin then [pkgs.darwin.apple_sdk.frameworks.Cocoa] else []);
+  executableSystemDepends = [
+    pkgs.blas
+    pkgs.liblapack
+    patched-gsl
+    sundials
+    pkgs.zlib
+  ] ++ (if pkgs.stdenv.isDarwin then [pkgs.darwin.apple_sdk.frameworks.Accelerate] else []);
+  license = stdenv.lib.licenses.bsd3;
+};
+
+haskellPackages = if compiler == "default"
+  then pkgs.haskellPackages
+  else pkgs.haskell.packages.${compiler};
+
+variant = if doBenchmark then pkgs.haskell.lib.doBenchmark else pkgs.lib.id;
+
+patched-hmatrix = pkgs.haskellPackages.hmatrix.overrideAttrs (oldAttrs: rec {
+  src = nixpkgs.fetchgit {
+    url = git://github.com/albertoruiz/hmatrix;
+    rev = "d83b17190029c11e3ab8b504e5cdc917f5863120";
+    sha256 = "11wr59wg21rky59j3kkd3ba6aqns9gkh0r1fnhwhn3fp7zfhanqn";
+  };
+  postUnpack = ''
+    sourceRoot=''${sourceRoot}/packages/base
+    echo Source root reset to ''${sourceRoot}
+  '';
+});
+
+patched-hmatrix-gsl = pkgs.haskellPackages.hmatrix-gsl.overrideAttrs (oldAttrs: rec {
+  src = nixpkgs.fetchgit {
+    url = git://github.com/albertoruiz/hmatrix;
+    rev = "d83b17190029c11e3ab8b504e5cdc917f5863120";
+    sha256 = "11wr59wg21rky59j3kkd3ba6aqns9gkh0r1fnhwhn3fp7zfhanqn";
+  };
+  postUnpack = ''
+    sourceRoot=''${sourceRoot}/packages/gsl
+    echo Source root reset to ''${sourceRoot}
+  '';
+});
+
+patched-gsl = pkgs.gsl.overrideAttrs (oldAttrs: rec {
+  src = nixpkgs.fetchgit {
+    url = git://github.com/idontgetoutmuch/gsl;
+    rev = "c2035977d65cd804169ff3370da6723cf879be75";
+    sha256 = "1fqp77gp9nl3av1z58cwg8fivik4rff394wgjzc76ayd04y0d1k7";
+    };
+  version = "2.5";
+  name = "gsl-${version}";
+  doCheck = false;
+  CFLAGS = "-DDEBUG";
+});
+
+drv = variant (haskellPackages.callPackage f {
+  hmatrix = patched-hmatrix;
+  hmatrix-gsl = patched-hmatrix-gsl;
+});
+
+in
+
+  if pkgs.lib.inNixShell then drv.env else drv
diff --git a/stack.yaml b/stack.yaml
index 1198e57..04bbfb6 100644
--- a/stack.yaml
+++ b/stack.yaml
@@ -13,5 +13,11 @@ packages:
 - packages/gsl/
 - packages/glpk/
 - packages/base/
-extra-deps: []
-resolver: lts-8.15
+- packages/sundials/
+extra-deps:
+- diagrams-rasterific-1.4
+- plots-0.1.0.2
+resolver: lts-10.9
+nix:
+  path: [nixpkgs=./nixpkgs.nix]
+  shell-file: shell.nix