packages feed

sde-solver (empty) → 0.1.0.0

raw patch · 10 files changed

+487/−0 lines, 10 filesdep +basedep +cerealdep +cereal-vectorsetup-changed

Dependencies added: base, cereal, cereal-vector, ghc-prim, haskell-mpi, mersenne-random-pure64, mtl, mwc-random, normaldistribution, parallel, vector

Files

+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2012, David Nilsson++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of David Nilsson nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ sde-solver.cabal view
@@ -0,0 +1,58 @@+name:                sde-solver++version:             0.1.0.0+synopsis:            Distributed SDE solver+description:         +  This package contains utilities for solving SDE instances in various ways.+  Basically an 'SDE' instance is solved using some 'SDESolver' working with some distribution mechanism.+  Results are gathered at the end point of the specified interval.+  .+  Included in the package are ways of doing distributed calculations over an MPI cluster,+  or optionally only using the local solver with built in parallelization.+  Two SDE instances have been implemented; geometric brownian motion and the Langevin equation,+  see the haddock documentation for an example.+  .+  The main interface is accessible through "Numeric.DSDE" which provides various way of solving generic problems.+  This module supports either local or distributed calculations in the IO monad and gathering the results as a distribution.+  Under the surface there is also a working pure implementation for monadic environments, using a pure Mersenne twister PRNG.+  .+  Internally there are several abstractions used when dealing with each component building up a solution.+  Given some 'SDE' and 'SDESolver' instances, it is also required to have some PRNG providing normally distributed numbers.+  This has been implemented over some specific monads and only results of type 'Double'.+  All of the internal components are written with polymorphism in mind, acting over some monad instance and generic result types in all cases.++homepage:            https://github.com/davnils/sde-solver++license:             BSD3+license-file:        LICENSE++author:              David Nilsson+maintainer:          nilsson.dd+code@gmail.com++category:            Math, Numerical+build-type:          Simple+cabal-version:       >=1.8++library+  exposed-modules:     Numeric.DSDE, Numeric.DSDE.Distribute, Numeric.DSDE.SDESolver, Numeric.DSDE.RNG, Numeric.DSDE.SDE.GeometricBrownian, Numeric.DSDE.SDE.Langevin, Numeric.DSDE.SDE+  +  build-depends:+    base >=4.5 && < 5,+    cereal >=0.3 && < 1,+    cereal-vector >=0.2 && < 1,+    ghc-prim >=0.2 && < 1,+    haskell-mpi >=1.2 && < 2,+    mersenne-random-pure64 >= 0.2.0.3 && < 1,+    mtl >=2.1 && < 3,+    mwc-random >=0.12 && < 1,+    normaldistribution >= 1.1.0.3 && < 2,+    parallel >= 3.2 && < 4,+    vector >=0.10 && < 1+  +  hs-source-dirs:      src+  ghc-options:         -O2+++source-repository head+  type: git+  location: git://github.com/davnils/sde-solver
+ src/Numeric/DSDE.hs view
@@ -0,0 +1,83 @@+{-# LANGUAGE ConstraintKinds #-}+--------------------------------------------------------------------+-- |+-- Module : Numeric.DSDE+--+-- Main entry point in the distributed SDE solver.+--+-- This module provides wrappers for the most common use cases, +-- and the user only needs to choose if to perform a distributed calculation+-- and supply the problem parameters.+--+-- Here follows an example of solving the Langevin equation over+-- the interval [0,2] with the result being written to a file on the main MPI node.+-- This produces a file with 1000 entries, corresponding to the final value+-- in all realisations.+--+-- @+--    let r = 1.0+--        sigma = 0.1+--        y_0 = 3.0+--        end = 0.2+--        stepSize = 0.01+--        runs = 1000+--    withSolver Milstein (distribute (Langevin r sigma) y_0 end stepSize runs) >>=+--      writeResult "out"+-- @++module Numeric.DSDE (+  -- * Abstract solver interface+  distribute,+  local,+  SDEResult(..),++  -- ** Supplied SDE Solvers +  Milstein(..),+  EulerMaruyama(..),++  -- * Helper functions+  withSolver, writeResult+) where++import qualified Data.Vector.Unboxed as V+import Numeric.DSDE.Distribute+import GHC.Conc+import Numeric.DSDE.RNG+import Numeric.DSDE.SDE+import Numeric.DSDE.SDESolver+import System.IO++-- | distribute: Perform a distributed MPI calculation over a set of hosts supplied through the MPI environment.+--   Each host is fully utilized by running a number of GHC threads equivalent to the number cores available.+--   Internally it uses the highly performant MWC PRNG over the IO monad.+--+--   local: Skips the MPI runtime and only runs locally on GHC threads.+distribute, local :: (SDE sde, SDESolver solver)+  => sde Double    -- ^ SDE Instance+  -> Double        -- ^ Initial y_0 value+  -> Double        -- ^ End of interval, evaluating over [0, end]+  -> Double        -- ^ Step size used by solver+  -> Int           -- ^ Number of realisations+  -> solver        -- ^ Solving method+  -> IO SDEResult  -- ^ Result ++distribute = evalWrapper [Distr MPI]+local = evalWrapper [] ++-- | Wrapper function providing access to the internal API.+evalWrapper :: (SDE sde, SDESolver solver) => [DistributeInstance IO] -> sde Double -> Double -> Double -> Double -> Int -> solver -> IO SDEResult+evalWrapper distr sde y_0 end step runs solver = do+  let rng = initialize+  cores <- getNumCapabilities+  evaluate (distr, Local cores) (sde, solver, rng, IP (End end) y_0 step runs)++-- | Use the specified solver with an IO action.+withSolver :: (SDESolver solver, Monad m) => solver -> (solver -> m a) ->  m a+withSolver = flip id++-- | Write a 'SDEResult' to the supplied filename, output as ASCII text.+writeResult :: FilePath -> SDEResult -> IO ()+writeResult file result = withFile file WriteMode $ output result+  where+  output (Scalar value count) h = hPutStr h $ unwords [show value, show count]+  output (Distribution samples) h = V.mapM_ (hPrint h) samples
+ src/Numeric/DSDE/Distribute.hs view
@@ -0,0 +1,176 @@+{-# LANGUAGE MultiParamTypeClasses, ConstraintKinds, BangPatterns,+             ExistentialQuantification, FlexibleInstances, DeriveGeneric,+             DefaultSignatures, ScopedTypeVariables, TupleSections,+             FlexibleContexts #-}++module Numeric.DSDE.Distribute where++import Control.Applicative ((<$>))+import Control.Concurrent+import Control.Concurrent.MVar+import Control.Monad.Identity hiding (mapM)+import Control.Monad.State+import Control.Parallel+import Control.Parallel.MPI.Simple +import Data.Foldable (fold, foldl')+import Data.Monoid+import Data.Serialize (Serialize(..))+import qualified Data.Vector.Unboxed as V+import Data.Vector.Serialize+import GHC.Generics (Generic)+import Numeric.DSDE.RNG+import Numeric.DSDE.SDE+import Numeric.DSDE.SDESolver+import Prelude hiding (sum, init, map)++-- | Wrapper used by all distributors supplied to the 'evaluate' function.+data DistributeInstance m = forall a. Distribute a m => Distr a++-- | Container describing the result produced by an SDE solution.+data SDEResult = Scalar !Double !Int             -- ^ Average of all end-point values, with the number of samples recorded.+               | Distribution !(V.Vector Double) -- ^ All end-point samples stored as an unboxed vector.+  deriving (Generic, Show)++-- | Monoid instance used when folding results from multiple sources.+instance Monoid SDEResult where+  mempty = Scalar 0 0+  Scalar a n         `mappend` Scalar b n'        =+    Scalar (elemSum / fromIntegral s) s+    where+    s = n + n'+    elemSum = a * fromIntegral n + b * fromIntegral n'++  Distribution v     `mappend` Distribution v'    = Distribution $ v V.++ v'+  Scalar _ _         `mappend` d@(Distribution _) = d+  d@(Distribution _) `mappend` Scalar _ _         = d++-- | Serialize instance used by MPI.+instance Serialize SDEResult++-- | Internal abstraction over the choice of specifying either the interval length or the number of steps.+data Accuracy = End Double+              | Steps Int+  deriving (Generic, Show)++-- | Serialize instance used by MPI.+instance Serialize Accuracy++-- | Set of parameters supplied to solve an SDE problem.+data InstanceParams = IP {+  accuracy :: !Accuracy,+  start :: !Double,+  deltat :: !Double,+  simulations :: !Int }+  deriving (Generic, Show)++-- | Serialize instance used by MPI.+instance Serialize InstanceParams++-- | MPI cluster distributor.+data MPI = MPI++-- | Local evaluation using GHC threads+data Local = Local Int++type SDEConstraint b c g m p = (SDE b, SDESolver c, Parameter p, RNGGen g m p)+type SDEInstance b c g m p = (b p, c, Maybe Int -> m g, InstanceParams)++-- | Type class indicating the ability to distribute data in some way.+--   Several distributors may be chained.+class Monad m => Distribute a m where+  -- | Inject an SDE instance into the context.+  inject :: SDEConstraint b c g m p =>+    a -> SDEInstance b c g m p -> m (SDEInstance b c g m p)+  -- | Remove an SDE result from the context.+  remove :: a -> SDEResult -> m SDEResult++-- | Type class indicating ability to solve an SDE problem.+class Execute a m p where+  execute :: SDEConstraint b c g m p =>  a -> SDEInstance b c g m p -> m SDEResult++-- | Type class indicating ability to perform a set of actions in an efficient way.+class Mappable m p where+  map' :: RNGGen g m p => (Int, Maybe Int -> m g) -> (g -> m b) -> [a] -> m [b] ++-- | Mappable instance for the IO monad. Work is divided using forkIO.+instance Mappable IO Double where+  map' (seed, rng) f l = do+    rand <- rng (Just seed)+    seeds <- mapM (\_ -> (Just . round <$> getRand rand) >>= rng) l+    mapM splitWork seeds >>= mapM takeMVar+    where+    splitWork rand = do+      var <- newEmptyMVar+      forkIO $ f rand >>= putMVar var+      return var++-- | Mappable instance for the pure State monad, uses 'par' annotations.+--   This does not perform well in general and needs to be optimized to+--   compete with the monadic IO instance.+instance RealFrac a => Mappable (State s) a where+  map' (seed, rng) f l = (go f seed l >>= sequence)+    where+    go f _ [] = return []+    go f s (_:t) = do+      worker <- rng (Just s)+      s' <- round <$> getRand worker+      rest <- go f s' t+      return $ f worker `par` f worker : rest++-- | Distribute instance over MPI which defines data transportation.+instance Distribute MPI IO where+  inject _ (sde, solver, rng, params) = do+    init+    size <- commSize commWorld+    rank <- commRank commWorld+    (sde, solver, rng,) <$> case rank of+      0 -> do+        let slaveSize = ceiling $ (fromIntegral $ simulations params :: Double) /+                                   fromIntegral size+        let slave = params { simulations = slaveSize}+        bcastSend commWorld 0 slave+        return slave++      _ -> bcastRecv commWorld 0++  remove _ localResult = do+    result <- commRank commWorld >>= retrieve+    finalize+    return result+    where+    retrieve 0 = do+      clusterResult <- gatherRecv commWorld 0 localResult+      return $ fold clusterResult+    retrieve _ =+      gatherSend commWorld 0 localResult >> return localResult++-- | Generic execute instances over any mappable monad 'm'.+instance Mappable m Double => Execute Local m Double where+  execute (Local cores) (!sde, !solver, !rng, params) = do+    seedRNG <- round <$> (rng Nothing >>= getRand)+    fold <$> map' (seedRNG, rng) runThread [1..cores]+    where+    perThread = ceiling $ (fromIntegral $ simulations params :: Double) /+                           fromIntegral cores :: Int+    steps = case accuracy params of+      End endTime -> floor $ endTime / deltat params+      Steps n -> n++    runThread rng = Distribution <$> thread rng+    thread rand = V.mapM (const $ threadEvaluation rand) $ V.replicate perThread (0.0 :: Double)+    threadEvaluation rand = foldM' (eval rand) (start params) [1..steps]+    eval rand w_i step = w_iplus1 solver sde rand (fromIntegral step * deltat params) w_i (deltat params)++-- | Evaluate the SDE using the supplied distributors and execution method.+evaluate :: (Monad m, SDEConstraint b c g m p, Execute e m p) =>+  ([DistributeInstance m], e) -> SDEInstance b c g m p -> m SDEResult+evaluate ([], method) input = execute method input+evaluate (Distr method : other, final) input =+  inject method input >>= evaluate (other, final) >>= remove method++-- | Monadic strict fold.+foldM' :: Monad m => (a -> b -> m a) -> a -> [b] -> m a+foldM' _ z [] = return z+foldM' func z (x:xs) = do+  z' <- func z x+  z' `seq` foldM' func z' xs
+ src/Numeric/DSDE/RNG.hs view
@@ -0,0 +1,51 @@+{-# LANGUAGE MultiParamTypeClasses, TypeSynonymInstances, FlexibleInstances,+             TypeFamilies, IncoherentInstances, ConstraintKinds,+             FunctionalDependencies, BangPatterns #-}++module Numeric.DSDE.RNG where++import Control.Applicative+import Control.Monad.Identity+import Control.Monad.ST+import Control.Monad.State+import Data.Maybe+import qualified Data.Random.Normal as NORMAL+import qualified Data.Vector.Unboxed as V+import qualified System.Random.Mersenne.Pure64 as MT+import Numeric.DSDE.SDE (Parameter)+import qualified System.Random.MWC as MWC+import qualified System.Random.MWC.Distributions as MD++-- | Default seed used by instances when the 'initialize' method is provided with 'Nothing'.+defaultSeed :: Int+defaultSeed = 0++-- | Typeclass describing a PRNG working in some monad 'm' and generating values of type 'p'.+class (Monad m, Functor m, Parameter p) => RNGGen g m p | g m -> p where+  -- | Generate a N(0,1) random number using the supplied generator.+  getRand :: g -> m p+  -- | Initialize a generator using some seed or defaulting to a constant.+  initialize :: Maybe Int -> m g++-- | 'RNGGen' instance for the ST-monadic MWC RNG.+instance RNGGen (MWC.Gen s) (ST s) Double where+  getRand = MD.normal 0 1+  initialize s = MWC.initialize . V.singleton . fromIntegral $ fromMaybe defaultSeed s++-- | 'RNGGen' instance for the IO-monadic MWC RNG.+instance (MWC.GenIO ~ d) => RNGGen d IO Double where+  {-# INLINE getRand #-}+  getRand = MD.normal 0 1+  {-# INLINE initialize  #-}+  initialize (Just n) = MWC.initialize . V.singleton . fromIntegral $ n+  initialize Nothing = MWC.withSystemRandom . MWC.asGenIO $ return++-- | 'RNGGen' instance for the pure monadic Mersenne Twister RNG, operating in monad 'State'.+instance RNGGen MT.PureMT (State MT.PureMT) Double where+  {-# INLINE getRand #-}+  getRand _ = do+    !(sample, gen') <- NORMAL.normal <$> get+    put gen'+    return sample+  {-# INLINE initialize  #-}+  initialize = return . MT.pureMT . fromIntegral . fromMaybe defaultSeed
+ src/Numeric/DSDE/SDE.hs view
@@ -0,0 +1,11 @@+{-# LANGUAGE ConstraintKinds #-}++module Numeric.DSDE.SDE where++-- | Type of result and intermediate variables used when solving an SDE problem.+type Parameter a = (Floating a, Fractional a, Num a)++-- | Describes an SDE problem. The partial derivate is used by advanced solvers+--   such as the 'Milstein' method.+class SDE a where+  f,g,partgoverparty :: Parameter p => a p -> p -> p -> p
+ src/Numeric/DSDE/SDE/GeometricBrownian.hs view
@@ -0,0 +1,16 @@+{-# Language BangPatterns #-}++module Numeric.DSDE.SDE.GeometricBrownian where++import Numeric.DSDE.SDE++-- | Geometric Brownian motion with rate r and diffusion sigma.+data GeometricBrownian p = GB !p !p++instance SDE GeometricBrownian where+  {-# SPECIALIZE INLINE f :: GeometricBrownian Double -> Double -> Double -> Double #-}+  f !(GB rate _) _ !w_i = rate *  w_i+  {-# SPECIALIZE INLINE g :: GeometricBrownian Double -> Double -> Double -> Double #-}+  g !(GB _ sigma) _ !w_i = sigma * w_i+  {-# SPECIALIZE INLINE partgoverparty :: GeometricBrownian Double -> Double -> Double -> Double #-}+  partgoverparty !(GB _ sigma) _ _ = sigma
+ src/Numeric/DSDE/SDE/Langevin.hs view
@@ -0,0 +1,16 @@+{-# Language BangPatterns #-}++module Numeric.DSDE.SDE.Langevin where++import Numeric.DSDE.SDE++-- | SDE instance of the Langevin equation with rate r and diffusion sigma.+data Langevin p = Langevin !p !p++instance SDE Langevin where+  {-# SPECIALIZE INLINE f :: Langevin Double -> Double -> Double -> Double #-}+  f !(Langevin r _) _ !w_i = -r * w_i+  {-# SPECIALIZE INLINE g :: Langevin Double -> Double -> Double -> Double #-}+  g !(Langevin _ sigma) _ _ = sigma+  {-# SPECIALIZE INLINE partgoverparty :: Langevin Double -> Double -> Double -> Double #-}+  partgoverparty _ _ _ = 0
+ src/Numeric/DSDE/SDESolver.hs view
@@ -0,0 +1,44 @@+{-# LANGUAGE MultiParamTypeClasses, TypeSynonymInstances, FlexibleInstances,+             BangPatterns, ConstraintKinds #-}++module Numeric.DSDE.SDESolver where++import Numeric.DSDE.SDE.GeometricBrownian+import Numeric.DSDE.RNG+import Numeric.DSDE.SDE+import qualified System.Random.MWC as M++-- | The Euler-Maruyama solving method. Order 1/2.+data EulerMaruyama = EulerMaruyama++-- | The Milstein solving method. Order 1.+data Milstein = Milstein++-- | Type class describing a method of solving SDE problems.+--    Defined by the next value produced in a solving sequence.+class SDESolver a where+  w_iplus1 :: (Monad m, SDE sde, RNGGen rng m p, Parameter p) =>+    a -> sde p -> rng -> p -> p -> p -> m p+  solverName :: a -> String++instance SDESolver EulerMaruyama where+  {-# INLINE w_iplus1 #-}+  {-# SPECIALIZE w_iplus1 :: EulerMaruyama -> GeometricBrownian Double -> M.GenIO -> Double -> Double -> Double -> IO Double #-}+  w_iplus1 _ !sde !rng !t_i !w_i !deltat = getRand rng >>= \rand -> return $+                            w_i+                            + f sde t_i w_i * deltat+                            + g sde t_i w_i * deltaB rand+    where deltaB r = sqrt deltat * r++  solverName _ = "Euler-Maruyama"++instance SDESolver Milstein where+  w_iplus1 _ !sde !rng !t_i !w_i !deltat = getRand rng >>= \rand -> return $+                            w_i+                            + f sde t_i w_i * deltat+                            + g' * deltaB rand+                            + g'/2 * partgoverparty sde t_i w_i  * (deltaB rand^^(2 :: Integer) - deltat)+    where+    deltaB r = sqrt deltat * r+    g' = g sde t_i w_i+  solverName _ = "Milstein"