Merge pull request #3 from FP-Modeling/refactor-readert

Refactor CT as a type alias for ReaderT

fact.cabal

@@ -54,6 +54,7 @@ library
     , criterion-measurement >=0.2.2.0
     , monadlist >=0.0.2
     , mtl >=2.3.1
+    , transformers >= 0.6.1.1
   default-language: Haskell2010
 
 executable fact-exe
@@ -70,5 +71,6 @@ executable fact-exe
     , criterion-measurement >=0.2.2.0
     , monadlist >=0.0.2
     , mtl >=2.3.1
+    , transformers >= 0.6.1.1    
     , fact    
   default-language: Haskell2010

src/Benchmarks.hs

@@ -2,6 +2,7 @@ module Benchmarks where
 
 import Examples.ChemicalReaction
 import Examples.Lorenz
+import Examples.Sine
 import Driver
 import CT
 import IO
@@ -18,7 +19,7 @@ import qualified Criterion.Measurement.Types as Criterion.Measurement.Measured
 
 perform :: IO [Double] -> IO (Double, Maybe Int64)
 perform test = do
-  (performance, _) <- measure (nfIO test) 10
+  (performance, _) <- measure (nfIO test) 1
   return (Criterion.Measurement.Measured.measTime performance, Criterion.Measurement.Measured.fromInt $ Criterion.Measurement.Measured.measAllocated performance)
 
 benchmarks :: IO ()

src/CT.hs

@@ -37,18 +37,17 @@
 -- Stability  : stable
 -- Tested with: GHC 8.10.7
 -- |
-
+{-# LANGUAGE FlexibleInstances #-}
 module CT 
        (CT(..),
         Parameters(..)) where
 
-import Control.Monad
-import Control.Monad.Fix
+import Control.Monad.Trans.Reader ( ReaderT )
 
-import Types
+import Types ( Iteration )
 
-import Solver
-import Simulation
+import Solver ( Solver )
+import Simulation ( Interval )
 
 -- | It defines the simulation time appended with additional information.
 data Parameters = Parameters { interval  :: Interval, -- ^ the simulation interval
@@ -57,83 +56,38 @@ data Parameters = Parameters { interval  :: Interval, -- ^ the simulation interv
                                iteration :: Iteration -- ^ the current iteration
                              } deriving (Eq, Show)
 
-newtype CT a = CT {apply :: Parameters -> IO a}
-
-instance Functor CT where
-  fmap f (CT da) = CT $ \ps -> fmap f (da ps)
-
-instance Applicative CT where
-  pure a = CT $ const (pure a)
-  (CT df) <*> (CT da) = CT $ \ps -> do f <- df ps
-                                       fmap f (da ps)
-
-appComposition :: CT (a -> b) -> CT a -> CT b
-appComposition (CT df) (CT da)
-  = CT $ \ps -> df ps >>= \f -> fmap f (da ps)
-
-instance Monad CT where
-  return = pure
-  (CT m) >>= k = CT $ \ps -> do a <- m ps
-                                k a `apply` ps
-
-instance MonadFix CT where
-  -- mfix :: (a -> m a) -> m a
-  mfix f = 
-    CT $ \ps -> mfix ((`apply` ps) . f)
-
-returnD :: a -> CT a
-returnD a = CT $ const (return a)
-
-bindD :: (a -> CT b ) -> CT a -> CT b
-bindD k (CT m) = 
-  CT $ \ps -> m ps >>= \a -> (\(CT m') -> m' ps) $ k a
-
-bindD' :: (a -> CT b ) -> CT a -> CT b
-bindD' k (CT m) = CT $ \ps -> do
-  a <- m ps
-  k a `apply` ps
-
-instance Eq (CT a) where
-  x == y = error "<< Can't compare dynamics >>" 
-
-instance Show (CT a) where
-  showsPrec _ x = showString "<< CT >>"
-
-unaryOP :: (a -> b) -> CT a -> CT b
-unaryOP = fmap
 
-binaryOP :: (a -> b -> c) -> CT a -> CT b -> CT c
-binaryOP func da db = fmap func da <*> db
+type CT a = ReaderT Parameters IO a
 
 instance (Num a) => Num (CT a) where
-  x + y = binaryOP (+) x y
-  x - y = binaryOP (-) x y
-  x * y = binaryOP (*) x y
-  negate = unaryOP negate
-  abs = unaryOP abs
-  signum = unaryOP signum
+  x + y = (+) <$> x <*> y
+  x - y = (-) <$> x <*> y
+  x * y = (*) <$> x <*> y
+  negate = fmap negate
+  abs = fmap abs
+  signum = fmap signum
   fromInteger i = return $ fromInteger i
 
 instance (Fractional a) => Fractional (CT a) where
-  x / y = binaryOP (/) x y
-  recip = unaryOP recip
+  x / y = (/) <$> x <*> y
+  recip = fmap recip
   fromRational t = return $ fromRational t
 
 instance (Floating a) => Floating (CT a) where
   pi = return pi
-  exp = unaryOP exp
-  log = unaryOP log
-  sqrt = unaryOP sqrt
-  x ** y = binaryOP (**) x y
-  sin = unaryOP sin
-  cos = unaryOP cos
-  tan = unaryOP tan
-  asin = unaryOP asin
-  acos = unaryOP acos
-  atan = unaryOP atan
-  sinh = unaryOP sinh
-  cosh = unaryOP cosh
-  tanh = unaryOP tanh
-  asinh = unaryOP asinh
-  acosh = unaryOP acosh
-  atanh = unaryOP atanh
+  exp = fmap exp
+  log = fmap log
+  sqrt = fmap sqrt
+  x ** y = (**) <$> x <*> y
+  sin = fmap sin
+  cos = fmap cos
+  tan = fmap tan
+  asin = fmap asin
+  acos = fmap acos
+  atan = fmap atan
+  sinh = fmap sinh
+  cosh = fmap cosh
+  tanh = fmap tanh
+  asinh = fmap asinh
+  acosh = fmap acosh
+  atanh = fmap atanh

src/Driver.hs

@@ -1,9 +1,12 @@
 module Driver where
 
 import CT
+    ( CT, Parameters(Parameters, solver, interval, time, iteration) )
 import Solver
+    ( Solver(stage, dt), Stage(SolverStage, Interpolate), iterToTime )
 import Simulation
-import Types
+    ( Interval(Interval), iterationHiBnd, iterationBnds )
+import Control.Monad.Trans.Reader ( ReaderT(runReaderT) )
 
 type Model a = CT (CT a)
 
@@ -12,43 +15,43 @@ epslon = 0.00001
 -- | Run the simulation and return the result in the last 
 -- time point using the specified simulation specs.
 runCTFinal :: Model a -> Double -> Solver -> IO a
-runCTFinal (CT m) t sl = 
-  do d <- m Parameters { interval = Interval 0 t,
-                         time = 0,
-                         iteration = 0,
-                         solver = sl { stage = SolverStage 0 }}
+runCTFinal m t sl = 
+  do d <- runReaderT m $ Parameters { interval = Interval 0 t,
+                                      time = 0,
+                                      iteration = 0,
+                                      solver = sl { stage = SolverStage 0 }}
      subRunCTFinal d t sl
 
 -- | Auxiliary functions to runCTFinal
 subRunCTFinal :: CT a -> Double -> Solver -> IO a
-subRunCTFinal (CT m) t sl =
-  do let iv = Interval 0 t
-         n = iterationHiBnd iv (dt sl)
-         disct = iterToTime iv sl n (SolverStage 0)
-         x = m Parameters { interval = iv,
-                            time = disct,
-                            iteration = n,
-                            solver = sl { stage = SolverStage 0 }}
-     if disct - t < epslon
-     then x
-     else m Parameters { interval = iv,
-                         time = t,
-                         iteration = n,
-                         solver = sl { stage = Interpolate }}
+subRunCTFinal m t sl = do
+  let iv = Interval 0 t
+      n = iterationHiBnd iv (dt sl)
+      disct = iterToTime iv sl n (SolverStage 0)
+      x = runReaderT m $ Parameters { interval = iv,
+                                      time = disct,
+                                      iteration = n,
+                                      solver = sl { stage = SolverStage 0 }}
+  if disct - t < epslon
+  then x
+  else runReaderT m $ Parameters { interval = iv,
+                                   time = t,
+                                   iteration = n,
+                                   solver = sl { stage = Interpolate }}
 
 -- | Run the simulation and return the results in all 
 -- integration time points using the specified simulation specs.
 runCT :: Model a -> Double -> Solver -> IO [a]
-runCT (CT m) t sl = do
-  d <- m Parameters { interval = Interval 0 t,
-                      time = 0,
-                      iteration = 0,
-                      solver = sl { stage = SolverStage 0}}
+runCT m t sl = do
+  d <- runReaderT m $ Parameters { interval = Interval 0 t,
+                                   time = 0,
+                                   iteration = 0,
+                                   solver = sl { stage = SolverStage 0}}
   sequence $ subRunCT d t sl
 
 -- | Auxiliary functions to runCT
 subRunCT :: CT a -> Double -> Solver -> [IO a]
-subRunCT (CT m) t sl = do
+subRunCT m t sl = do
   let iv = Interval 0 t
       (nl, nu) = iterationBnds iv (dt sl)
       parameterize n =
@@ -65,8 +68,8 @@ subRunCT (CT m) t sl = do
                    iteration = nu,
                    solver = sl {stage = Interpolate}}
       endTime = iterToTime iv sl nu (SolverStage 0)
-      values = map (m . parameterize) [nl .. nu]
+      values = map (runReaderT m . parameterize) [nl .. nu]
   if endTime - t < epslon
   then values
-  else init values ++ [m ps]
+  else init values ++ [runReaderT m ps]
 

src/Examples/Sine.hs

@@ -1,4 +1,3 @@
-{-# LANGUAGE RecordWildCards #-}
 module Examples.Sine where
 
 import Driver
@@ -11,8 +10,8 @@ import Data.List
 import Simulation
 
 sineSolv = Solver { dt = 0.01,
-                method = RungeKutta4,
-                stage = SolverStage 0 }
+                    method = RungeKutta4,
+                    stage = SolverStage 0 }
 
 sineModel :: Model [Double]
 sineModel =

src/Examples/Test.hs

@@ -61,25 +61,25 @@ predicate initialCondition =
   if initialCondition >= 20
   then example2 else example1
 
-demux :: Predicate Double Result -> HybridModel Result
-demux predicate (initialCondition, _) p = do
-  let m = predicate initialCondition
-  model <- m (pure initialCondition) `apply` p
-  head <$> model `apply` p
+-- demux :: Predicate Double Result -> HybridModel Result
+-- demux predicate (initialCondition, _) p = do
+--   let m = predicate initialCondition
+--   model <- m (pure initialCondition) `apply` p
+--   head <$> model `apply` p
 
-hybrid :: (MonadPlus p, Monad m) => (a -> Parameters -> m a) -> a -> Double -> Solver -> m (p a)
-hybrid f z t sl =
-  do let iv = Interval 0 t
-         (nl, nu) = iterationBnds iv (dt sl)
-         parameterise n = Parameters { interval = iv,
-                                       time = iterToTime iv sl n (SolverStage 0),
-                                       iteration = 1,
-                                       solver = sl { stage = SolverStage 0 }}
-         ps = map parameterise [nl..nu]
-     scanM f z ps
+-- hybrid :: (MonadPlus p, Monad m) => (a -> Parameters -> m a) -> a -> Double -> Solver -> m (p a)
+-- hybrid f z t sl =
+--   do let iv = Interval 0 t
+--          (nl, nu) = iterationBnds iv (dt sl)
+--          parameterise n = Parameters { interval = iv,
+--                                        time = iterToTime iv sl n (SolverStage 0),
+--                                        iteration = 1,
+--                                        solver = sl { stage = SolverStage 0 }}
+--          ps = map parameterise [nl..nu]
+--      scanM f z ps
 
-test = do
-  t <- hybrid (demux predicate) (1, "initial") 40 sineSolv2
-  case t of
-    [] -> fail "Something went wrong during hybrid simulation"
-    list -> print list
+-- test = do
+--   t <- hybrid (demux predicate) (1, "initial") 40 sineSolv2
+--   case t of
+--     [] -> fail "Something went wrong during hybrid simulation"
+--     list -> print list

src/IO.hs

@@ -2,8 +2,8 @@ module IO where
 
 import System.IO
     ( Handle, hClose, openFile, hPutStrLn, IOMode(WriteMode) )
-import Simulation
-import Solver
+import Simulation ( Interval(stopTime, startTime) )
+import Solver ( Solver(dt) )
 
 addTime :: IO [[a]] -> Interval -> Solver -> IO [(Double, [a])]
 addTime answers intv solver = fmap (zip input) answers