From fb80f07072991ca1bddae0a19a266abbfb0a8341 Mon Sep 17 00:00:00 2001 From: Dominic Steinitz Date: Wed, 28 Mar 2018 11:20:23 +0100 Subject: Using hsc, semi-auto-generating documentation --- packages/sundials/hmatrix-sundials.cabal | 18 ++- packages/sundials/src/Bar.hsc | 6 +- packages/sundials/src/Main.hs | 32 +++++ .../sundials/src/Numeric/Sundials/ARKode/ODE.hs | 142 +++++++++++++++++++-- packages/sundials/stack.yaml | 66 ---------- 5 files changed, 180 insertions(+), 84 deletions(-) delete mode 100644 packages/sundials/stack.yaml (limited to 'packages') 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 -- cgit v1.2.3