packages feed

goal-simulation (empty) → 0.1

raw patch · 22 files changed

+2067/−0 lines, 22 filesdep +basedep +cairodep +clocksetup-changed

Dependencies added: base, cairo, clock, directory, goal-core, goal-geometry, goal-probability, goal-simulation, gtk, hmatrix, machines, mtl, vector

Files

+ Goal/Simulation.hs view
@@ -0,0 +1,30 @@+{-| A general purpose library for simulating differential processes, of a deterministic or+    stochastic nature. -}++module Goal.Simulation+    ( -- * Exports+      module Goal.Simulation.Mealy+    , module Goal.Simulation.Flow+    , module Goal.Simulation.Filter+    , module Goal.Simulation.Filter.Flow+    , module Goal.Simulation.Optimization+    , module Goal.Simulation.Chain+    , module Goal.Simulation.Plot+    , module Goal.Simulation.Physics.Configuration+    , module Goal.Simulation.Physics.Models.Pendulum+    ) where++--- Imports ---+++-- Goal --++import Goal.Simulation.Mealy+import Goal.Simulation.Optimization+import Goal.Simulation.Flow+import Goal.Simulation.Filter+import Goal.Simulation.Filter.Flow+import Goal.Simulation.Chain+import Goal.Simulation.Plot+import Goal.Simulation.Physics.Configuration+import Goal.Simulation.Physics.Models.Pendulum
+ Goal/Simulation/Chain.hs view
@@ -0,0 +1,70 @@+-- | A general purpose library for simulating differential processes, of a deterministic or+-- stochastic nature.++module Goal.Simulation.Chain (+    -- * Chains+      Chain+    , streamChain+    -- ** Markov Chains+    , MarkovTensor+    , markovTensor+    , markovChain+    -- ** Harris Chains+    , harrisChain+    ) where++--- Imports ---+++-- Goal --++import Goal.Simulation.Mealy++import Goal.Geometry+import Goal.Probability+++--- Stochastic Processes ---+++-- | A 'Chain' is a discrete time 'Flow'.+type Chain s = Mealy () s++streamChain :: Chain k -> [k]+-- | A convenience function for streaming 'MarkovChain's.+streamChain mchn = stream mchn $ repeat ()++-- Finite --++type MarkovTensor s = Tensor (CurvedCategorical s) (CurvedCategorical s)++markovTensor :: s -> MarkovTensor s+markovTensor s = Tensor (CurvedCategorical s) (CurvedCategorical s)++markovChain :: Discrete s+    => (Function Standard Standard :#: MarkovTensor s) -- ^ The stochastic matrix+    -> Element s -- ^ The initial state+    -> RandST r (Chain (Element s)) -- ^ An embedded markov chain+-- | Constructs a markov chain process.+markovChain smtx =+    accumulateRandomFunction (markovChainAccumulator smtx)++markovChainAccumulator smtx _ k = do+    let ks = elements . sampleSpace . domain $ manifold smtx+        sk = fromList (domain $ manifold smtx) [ if k' == k then 1 else 0 | k' <- ks ]+    k' <- generate $ smtx >.> sk+    return (k',k')++-- Continuous --++harrisChain :: (Manifold m, Transition c Standard m, Generative Standard m)+    => (Sample m -> c :#: m)+    -> Sample m+    -> RandST s' (Chain (Sample m))+harrisChain mdl =+    accumulateRandomFunction (harrisChainAccumulator mdl)++harrisChainAccumulator mdl () s = do+    let cp' = mdl s+    s' <- standardGenerate cp'+    return (s',s')
+ Goal/Simulation/Filter.hs view
@@ -0,0 +1,77 @@+-- | A general purpose library for simulating differential processes, of a deterministic or+-- stochastic nature.++module Goal.Simulation.Filter+    (+    -- * Learning+      beliefBackpropagation+    -- * Statistics+    , beliefNegativeLogLikelihoods+    -- * Structural+    , decomposeBeliefTransition+    ) where+++--- Imports ---+++-- Goal --++import Goal.Core++import Goal.Geometry+import Goal.Probability+++--- Learning ---++beliefBackpropagation+    :: ( Manifold x, ExponentialFamily m, ExponentialFamily y, ExponentialFamily n+       , Generative Natural m, Generative Natural n+       , Riemannian Natural m, Riemannian Natural y, Legendre Natural m, Legendre Natural y)+    => [NaturalFunction :#: Harmonium m n]+    -> Int+    -> Int+    -> (Function Mixture Mixture :#: NeuralNetwork m y m)+    -> [(c, Mixture, Mixture) :#: (x, m, m)]+    -> Rand (ST s) (Differentials :#: Tangent (Function Mixture Mixture) (NeuralNetwork m y m))+-- | NB: Assumes null latent state bias on the Harmoniums at the moment.+beliefBackpropagation trnss blkn cdn nnp xnzs = do+    let zns' = decomposeBeliefTransition <$> zip xnzs (tail xnzs)+        (zs,ns') = unzip zns'+        (eys,ys,ez0s,z0s) = feedForward nnp zs+        trnss' = harmoniumTranspose <$> zipWith modulateHarmoniumBelief z0s trnss+    gbss <- zipWithM (bulkGibbsSampling0 cdn) trnss' $ replicate blkn <$> ns'+    let project trns' z0 n' = potentialMapping $ trns' >.> (z0 <+> n')+        oNss = snd . unzip . last <$> gbss+        (_,mtx,_) = splitHarmonium $ head trnss+        errs = mtx >$> [ meanPoint [ project trns' z0 oN | oN <- oNs ] <-> project trns' z0 n'+                      | (trns',z0,oNs,n') <- zip4 trnss' z0s oNss ns' ]+    return $ feedBackward nnp zs eys ys ez0s errs++beliefNegativeLogLikelihoods+    :: (Manifold x, Manifold m, AbsolutelyContinuous Natural n, ExponentialFamily n, Sample n ~ [Double])+    => NaturalFunction :#: Harmonium m n+    -> (c, Mixture, Mixture) :#: (x, m, m)+    -> (Double, Double)+beliefNegativeLogLikelihoods trns xnz =+    let (x,n,z) = splitTriple xnz+        zx = harmoniumTranspose trns >.> z+        nx = harmoniumTranspose trns >.> n+     in (negate . log . density nx $ listCoordinates x, negate . log . density zx $ listCoordinates x)++++--- Internal ---+++decomposeBeliefTransition+    :: (Manifold m, Manifold n)+    => ( (c, Mixture, Mixture) :#: (m, n, n)+       , (c, Mixture, Mixture) :#: (m, n, n) )+    -> (Mixture :#: n, Mixture :#: n)+decomposeBeliefTransition (xnz,xnz') =+    let (_,_,z) = splitTriple xnz+        (_,n',_) = splitTriple xnz'+     in (z, n')+
+ Goal/Simulation/Filter/Flow.hs view
@@ -0,0 +1,97 @@+-- | A general purpose library for simulating differential processes, of a deterministic or+-- stochastic nature.++module Goal.Simulation.Filter.Flow+    (+    -- * Belief Processes+      coxProcess+    , filteringProcess+    , transducerProcess+    -- ** Statistics+    , locationBeliefField+    ) where+++--- Imports ---+++-- Goal --++import Goal.Core++import Goal.Simulation.Mealy++import Goal.Geometry+import Goal.Probability+import Goal.Simulation.Flow+++--- Processes ---+++coxProcess :: (ExponentialFamily m, Generative Natural m, Sample m ~ [Double])+    => (Double -> Double) -- ^ Gain over time+    -> (NaturalFunction :#: Harmonium (Replicated Poisson) m)+    -> RandST s (Mealy (Double, [Double]) [Int])+-- | Constructs a inhomogeneous poisson process. Since we have yet to+-- develop tools for describing distributions on manifolds directly, the+-- dynamic response works simply on sample spaces.+coxProcess gnf trns =+    accumulateRandomFunction (coxAccumulator gnf trns) 0++coxAccumulator gnf trns (t',cs) t = do+    let dt = t' - t+        trns' = modulateTransducerGain (dt*gnf t') trns+    n' <- generate $ conditionalLatentDistribution trns' cs+    return (n',t')++transducerProcess :: Manifold m+    => [Double]+    -> (Double -> Double)+    -> NaturalFunction :#: Harmonium (Replicated Poisson) m+    -> [NaturalFunction :#: Harmonium (Replicated Poisson) m]+-- | Right now assumes that dt is constant. Watch out for that.+transducerProcess ts@(t0:t1:_) gnf trns =+    let dt = t1 - t0+     in  [ modulateTransducerGain (dt * gnf t) trns | t <- ts ]+transducerProcess _ _ _ = error "Try Harder"+++filteringProcess+    :: (ExponentialFamily m, Apply Mixture Mixture f, Domain f ~ m, Codomain f ~ m, Manifold x)+    => Flow c x -- ^ Underlying Flow+    -> Mealy (Double, [Double]) (Sample m) -- ^ Emission Process+    -> (Function Mixture Mixture :#: f) -- ^ Belief Dynamics+    -> (Mixture :#: m) -- ^ Initial Belief+    -> Flow (c, Mixture, Mixture) (x, m, m) -- ^ The belief process+filteringProcess flow cox nnp z0 =+    accumulateMealy z0 $ proc (t',z) -> do+        cm' <- flow -< t'+        n' <- cox -< (t', listCoordinates cm')+        let mn' = sufficientStatistic (manifold z) n'+            z' = nnp >.> z <+> mn'+            p' = joinTriple cm' mn' z'+        returnA -< (p', z')+++--- Plot Oriented ---++locationBeliefField+    :: ( ExponentialFamily m, ExponentialFamily n, Generative Natural n+       , Apply Mixture Mixture f, Domain f ~ m, Codomain f ~ m+       , SampleSpace n ~ Continuum, Transition Natural Standard n )+    => NaturalFunction :#: Harmonium m (Replicated n)+    -> Double+    -> Double+    -> Int+    -> Int+    -> [Double]+    -> Function Mixture Mixture :#: f+    -> (Double,Double)+    -> (Double,Double)+locationBeliefField trns dt scl ix iy xs nnp (x,y) =+    let xs' = take ix xs ++ x : drop (ix+1) xs+        xs'' = take iy xs' ++ y : drop (iy+1) xs'+        sp = chart Standard . transition . (harmoniumTranspose trns >.>)+            . (nnp >.>) . potentialMapping $ conditionalLatentDistribution trns xs''+     in (scl/dt * (coordinate (2*ix) sp - (xs'' !! ix)), scl/dt * (coordinate (2*iy) sp - (xs'' !! iy)))
+ Goal/Simulation/Flow.hs view
@@ -0,0 +1,194 @@+-- | A general purpose library for simulating differential processes, of a deterministic or+-- stochastic nature.++module Goal.Simulation.Flow (+    -- * Flows+      Flow+    -- ** Deterministic+    , autonomousODE+    , nonAutonomousODE+    , lagrangianFlow+    -- ** Stochastic+    , itoProcess+    , langevinStep+    , langevinFlow+    -- * Integral Curves+    , stepEuler+    , stepRK4+    , stepEuler'+    , stepRK4'+    ) where++--- Imports ---+++-- Goal --++import Goal.Core++import Goal.Simulation.Mealy++import Goal.Geometry+import Goal.Simulation.Physics.Configuration+import Goal.Probability++-- Qualified --++import qualified Numeric.LinearAlgebra.HMatrix as M+++--- Deterministic ---+++-- | A 'Flow' can be used to generate a path of 'Element's through a particular+-- space 's' in time ('Double').+type Flow c m = Mealy Double (c :#: m)++autonomousODE :: Manifold m+    => (Coordinates -> Coordinates) -- ^ Differential Equation+    -> (c :#: m) -- ^ Initial State+    -> Flow c m -- ^ Differential Process+-- | Creates a process out of the ordinary differential equation described by the given+-- explicit function.+autonomousODE f' p0 =+    accumulateFunction accumulator (0,p0)+      where accumulator t' (t,p)+                | t' == t = (p,(t,p))+                | otherwise =+                      let dt = t' - t+                          p' = p <+> fromCoordinates (manifold p) (stepRK4 f' dt $ coordinates p)+                       in (p',(t',p'))++nonAutonomousODE :: Manifold m+    => (Double -> Coordinates -> Coordinates) -- ^ Differential Equation+    -> Double -- ^ Initial Time+    -> (c :#: m) -- ^ Initial State+    -> Flow c m -- ^ Differential Process+-- | Creates a process out of the ordinary differential equation described by the given+-- explicit function.+nonAutonomousODE f' t0 p0 =+    accumulateFunction accumulator (t0,p0)+      where accumulator t' (t,p)+                | t' == t = (p,(t,p))+                | otherwise =+                      let dt = t' - t+                          p' = p <+> fromCoordinates (manifold p) (stepRK4' f' t dt $ coordinates p)+                       in (p',(t',p'))++-- Stochastic --++itoProcess :: Manifold m+    => (Double -> Coordinates -> Coordinates) -- ^ The drift function+    -> (Double -> Coordinates -> M.Matrix Double) -- ^ The diffusion function+    -> Double -- ^ The initial time+    -> (c :#: m) -- ^ The initial state+    -> RandST s (Flow c m) -- ^ The Ito flow+-- | Constructs an ito process.+itoProcess mu sgma t0 p0 =+    accumulateRandomFunction (itoAccumulator mu sgma) (t0,p0)++itoAccumulator mu sgma t' (t,p)+  | t' == t = return (p,(t,p))+  | otherwise = do+        let dt = t' - t+            x = coordinates p+            x' = x + stepRK4' mu t dt x+        x'' <- generate $ muSigmaToMultivariateNormal x' (M.scale (sqrt dt) $ sgma t x)+        let p' = fromCoordinates (manifold p) x''+        return  (p',(t',p'))++-- Mechanical --++lagrangianFlow :: (Riemannian Generalized m, ForceField f m)+    => f+    -> (Partials :#: PhaseSpace m)+    -> Flow Partials (PhaseSpace m)+lagrangianFlow f p0 =+    accumulateFunction accumulator (0,p0)+      where accumulator t' (t,qdq)+                | t' == t = (qdq,(t,qdq))+                | otherwise =+                      let dt = t' - t+                          ddq xs = coordinates . vectorField f $ fromCoordinates (manifold qdq) xs+                          qdq' = qdq <+> fromCoordinates (manifold qdq) (stepRK4 ddq dt $ coordinates qdq)+                       in (qdq', (t',qdq'))++langevinStep :: (Riemannian Generalized m, ForceField f m)+    => Double+    -> f+    -> ( Partials :#: PhaseSpace m+       -> Function Differentials Differentials :#: Tensor (GeneralizedAcceleration m) (GeneralizedAcceleration m) )+    -> (Partials :#: PhaseSpace m)+    -> RandST s (Partials :#: PhaseSpace m)+langevinStep dt f sgma qdq = do+    let flx p = matrixSquareRoot $ sgma p+        dq = bundleToTangent qdq+    nrms <- replicateM (dimension $ manifold dq) . generate . chart Standard $ fromList Normal [0,1]+    let lng p = sqrt dt /> (flx p >.> fromList (Tangent $ bundleToTangent p) nrms)+        ddq xs = coordinates . vectorField (f, lng) $ fromCoordinates (manifold qdq) xs+        qdq' = qdq <+> fromCoordinates (manifold qdq) (stepRK4 ddq dt $ coordinates qdq)+    return qdq'++langevinFlow :: (Riemannian Generalized m, ForceField f m)+    => f+    -> ( Partials :#: PhaseSpace m+           -> Function Differentials Differentials :#: Tensor (GeneralizedAcceleration m) (GeneralizedAcceleration m) )+    -> (Partials :#: PhaseSpace m)+    -> RandST s (Flow Partials (PhaseSpace m))+langevinFlow f sgma p0 =+    accumulateRandomFunction (langevinAccumulator f sgma) (0,p0)++langevinAccumulator f sgma t' (t,qdq)+    | t' == t = return (qdq,(t,qdq))+    | otherwise = do+          let dt = t' - t+          qdq' <- langevinStep dt f sgma qdq+          return (qdq', (t',qdq'))+++--- Integral Curves ---+++stepEuler+    :: (Coordinates -> Coordinates) -- ^ The derivative of the function to simulate+    -> Double -- ^ The time step 'dt'+    -> Coordinates -- ^ The state of the system at the current time+    -> Coordinates -- ^ The difference to the state of the system at time 't' + 'dt'+-- | Returns the difference from an Euler step.+stepEuler f' dt x = realToFrac dt * f' x++stepRK4+    :: (Coordinates -> Coordinates) -- ^ The derivative of the function to simulate+    -> Double -- ^ The time step 'dt'+    -> Coordinates -- ^ The state of the system at the current time+    -> Coordinates -- ^ The difference to the state of the system at time 't' + 'dt'+-- | Returns the difference from an RK4 step.+stepRK4 f' dt x =+    let k1 = realToFrac dt * f' x+        k2 = realToFrac dt * f' (x + k1 / 2)+        k3 = realToFrac dt * f' (x + k2 / 2)+        k4 = realToFrac dt * f' (x + k3)+    in (k1 + 2 * k2 + 2 * k3 + k4) / 6++stepEuler'+    :: (Double -> Coordinates -> Coordinates) -- ^ The derivative of the function to simulate+    -> Double -- ^ The current time 't'+    -> Double -- ^ The time step 'dt'+    -> Coordinates -- ^ The state of the system at the current time+    -> Coordinates -- ^ The difference to the state of the system at time 't' + 'dt'+-- | Time inhomogenous version.+stepEuler' f' t dt x = realToFrac dt * f' t x++stepRK4'+    :: (Double -> Coordinates -> Coordinates) -- ^ The derivative of the function to simulate+    -> Double -- ^ The current time 't'+    -> Double -- ^ The time step 'dt'+    -> Coordinates -- ^ The state of the system at the current time+    -> Coordinates -- ^ The difference to the state of the system at time 't' + 'dt'+-- | Time inhomogenous version.+stepRK4' f' t dt x =+    let k1 = realToFrac dt * f' t x+        k2 = realToFrac dt * f' (t + dt / 2) (x + k1 / 2)+        k3 = realToFrac dt * f' (t + dt / 2) (x + k2 / 2)+        k4 = realToFrac dt * f' (t + dt) (x + k3)+    in (k1 + 2 * k2 + 2 * k3 + k4) / 6
+ Goal/Simulation/Mealy.hs view
@@ -0,0 +1,119 @@+{-| A general purpose library for simulating differential processes, of a deterministic or+    stochastic nature. -}++module Goal.Simulation.Mealy+    ( -- * Exports+      module Data.Machine+    -- * Accumulation+    , accumulateFunction+    , accumulateFunction'+    , accumulateMealy+    , accumulateMealy'+    , accumulateRandomFunction+    , accumulateRandomFunction'+    , accumulateRandomFunction0+      -- * Execution+    , stream+    , streamM+    , streamM_+    ) where++--- Imports ---+++-- Goal --++import Goal.Probability++-- Reexporting --++import Data.Machine++import qualified Control.Monad.ST as ST+++--- Mealys ---+++accumulateFunction :: (a -> acc -> (b,acc)) -> acc -> Mealy a b+-- | accumulateFunction takes a function from a value and an accumulator (e.g. just a sum+-- value or an evolving set of parameters for some model) to a value and an accumulator.+-- The accumulator is then looped back into the function, returning a Mealy from a to+-- b, which updates the accumulator every time step.+accumulateFunction f acc = Mealy $ \a ->+    let (b,acc') = f a acc+     in (b,accumulateFunction f acc')++accumulateFunction' :: (a -> acc -> (b,acc)) -> acc -> Mealy a (b,acc)+-- | accumulateFunction' acts like accumulateFunction but the Mealy automata will+-- continue to return the accumulator as it generates it.+accumulateFunction' f =+    accumulateFunction f'+    where f' a acc =+              let (b,acc') = f a acc+               in ((b,acc'),acc')++accumulateRandomFunction :: (a -> acc -> forall s . RandST s (b,acc)) -> acc -> RandST s' (Mealy a b)+-- | accumulateRandomFunction is analogous to accumulateFunction, but takes as an+-- argument a function which returns a random variable.+accumulateRandomFunction rf acc0 = do+    rf' <-  accumulateRandomFunction0 (uncurry rf)+    return $ accumulateMealy acc0 rf'++accumulateRandomFunction' :: (a -> acc -> forall s . RandST s (b,acc)) -> acc -> RandST s' (Mealy a (b,acc))+-- | accumulateRandomFunction' is analogous to accumulateFunction', but takes as an+-- argument a function which returns a random variable.+accumulateRandomFunction' rf acc0 = do+    rf' <- accumulateRandomFunction0 (uncurry rf)+    return $ accumulateMealy' acc0 rf'++accumulateRandomFunction0 :: (a -> forall s . RandST s b) -> RandST s' (Mealy a b)+-- | accumulateRandomFunction' Mealifies stateless random functions.+accumulateRandomFunction0 rf = do+    sd <- seed+    return $ accumulateFunction f sd+    where f a sd = ST.runST $ do+              gn <- restore sd+              b <- runRand (rf a) gn+              sd' <- save gn+              return (b,sd')++accumulateMealy :: acc -> Mealy (a,acc) (b,acc) -> Mealy a b+-- | accumulateMealy takes a Mealy with an accumulating parameter and loops it.+accumulateMealy acc0 mly0 =+    accumulateFunction f (acc0,mly0)+    where f a (acc,Mealy cf) =+              let ((b,acc'),mly') = cf (a,acc)+               in (b,(acc',mly'))++accumulateMealy' :: acc -> Mealy (a,acc) (b,acc) -> Mealy a (b,acc)+-- | accumulateMealy except with a returned accumulator.+accumulateMealy' acc0 mly0 =+    accumulateFunction f (acc0,mly0)+    where f a (acc,Mealy cf) =+              let ((b,acc'),mly') = cf (a,acc)+               in ((b,acc'),(acc',mly'))++--- Execution ---++{-+parallelizeMealys :: [Mealy a b] -> Mealy [a] [b]+{-| Turns a list of circuits into a circuit over lists, bound by the power of parMap rseq. -}+parallelizeMealys crcs = Mealy $ \as ->+    let (bs,crcs') = unzip $ parZip crcs as+    in  (bs, parallelizeMealys crcs')+    where parZip [] [] = []+          parZip (crc:crcs) (a:as) =+              let (b,crc') = runMealy crc a+              in  b `par` (b,crc') : parZip crcs as+          parZip _ _ = error "Parallel circuit does not match size of input"+          -}++stream :: Mealy a b -> [a] -> [b]+stream mly as = run . supply as . auto $ mly++streamM :: Monad m => Mealy a b -> (b -> m c) -> [a] -> m [c]+streamM mly fM as = runT . supply as $ auto mly ~> autoM fM++streamM_ :: Monad m => Mealy a b -> (b -> m c) -> [a] -> m ()+streamM_ mly fM as = runT_ . supply as $ auto mly ~> autoM fM
+ Goal/Simulation/Optimization.hs view
@@ -0,0 +1,200 @@+-- | The Map module provides tools for developing function space 'Manifold's.+-- A map is a 'Manifold' where the 'Point's of the Manifold represent+-- parametric functions between 'Manifold's. The defining feature of 'Map's is+-- that they have a particular 'Domain' and 'Codomain', which themselves are+-- 'Manifold's.++module Goal.Simulation.Optimization (+    -- * Mean Squared Error+      meanSquaredError+    -- * Cauchy Sequences+    , cauchyLimit+    , cauchySequence+    -- * Gradient Pursuit+    , stochasticGradientDescent+    , stochasticGradientAscent+    , stochasticVanillaGradientDescent+    , stochasticVanillaGradientAscent+    , boundedStochasticVanillaGradientDescent+    , boundedStochasticVanillaGradientAscent+    -- * Least Squares+    , designMatrix+    , leastSquares+    , leastSquares0+    -- ** Newton+    , newtonStep+    , newtonSequence+    -- ** Gauss Newton+    , gaussNewtonStep+    ) where++--- Imports ---++import Prelude hiding (map,minimum,maximum)++-- Goal --++import Goal.Core+import Goal.Geometry+import Goal.Probability++import Goal.Simulation.Mealy+++--- Stochastic Pursuit ---+++type StochasticPursuit x c m = Mealy x (c :#: m)+++--- Gradient Descent ---++stochasticGradientAscent :: (Riemannian c m, Manifold m)+    => Double -- ^ Step size+    -> (c :#: m -> x -> forall s . RandST s (Differentials :#: Tangent c m)) -- ^ Gradient calculator+    -> (c :#: m) -- ^ The initial point+    -> RandST r (StochasticPursuit x c m) -- ^ The gradient ascent+stochasticGradientAscent eps0 f0' = accumulateRandomFunction (accumulator eps0 f0')+    where accumulator eps f' x cm = do+              dcm <- f' cm x+              let cm' = gradientStep eps $ sharp dcm+              return (cm',cm')++boundedStochasticVanillaGradientAscent :: Manifold m+    => Double -- ^ Step size+    -> Double -- ^ Gradient Rejection Bound+    -> (c :#: m -> x -> forall s . RandST s (Differentials :#: Tangent c m)) -- ^ Gradient calculator+    -> (c :#: m) -- ^ The initial point+    -> RandST r (StochasticPursuit x c m) -- ^ The gradient ascent+boundedStochasticVanillaGradientAscent eps0 bnd0 f0' = accumulateRandomFunction (accumulator eps0 bnd0 f0')+    where accumulator eps bnd f' x cm = do+              dcm <- f' cm x+              let cm' = if maximum (abs <$> listCoordinates dcm) < bnd+                            then gradientStep eps $ breakChart dcm+                            else trace "Ping!" cm+              return (cm',cm')++stochasticVanillaGradientAscent :: Manifold m+    => Double -- ^ Step size+    -> (c :#: m -> x -> forall s . RandST s (Differentials :#: Tangent c m)) -- ^ Gradient calculator+    -> (c :#: m) -- ^ The initial point+    -> RandST r (StochasticPursuit x c m) -- ^ The gradient ascent+stochasticVanillaGradientAscent eps0 f0' = accumulateRandomFunction (accumulator eps0 f0')+    where accumulator eps f' x cm = do+              dcm <- f' cm x+              let cm' = gradientStep eps $ breakChart dcm+              return (cm',cm')++stochasticGradientDescent :: (Riemannian c m, Manifold m)+    => Double -- ^ Step size+    -> (c :#: m -> x -> forall s . RandST s (Differentials :#: Tangent c m)) -- ^ Gradient calculator+    -> (c :#: m) -- ^ The initial point+    -> RandST r (StochasticPursuit x c m) -- ^ The gradient ascent+stochasticGradientDescent eps = stochasticGradientAscent (-eps)++boundedStochasticVanillaGradientDescent :: Manifold m+    => Double -- ^ Step size+    -> Double -- ^ Gradient Rejection Bound+    -> (c :#: m -> x -> forall s . RandST s (Differentials :#: Tangent c m)) -- ^ Gradient calculator+    -> (c :#: m) -- ^ The initial point+    -> RandST r (StochasticPursuit x c m) -- ^ The gradient ascent+boundedStochasticVanillaGradientDescent eps = boundedStochasticVanillaGradientAscent (-eps)++stochasticVanillaGradientDescent :: Manifold m+    => Double -- ^ Step size+    -> (c :#: m -> x -> forall s . RandST s (Differentials :#: Tangent c m)) -- ^ Gradient calculator+    -> (c :#: m) -- ^ The initial point+    -> RandST r (StochasticPursuit x c m) -- ^ The gradient ascent+stochasticVanillaGradientDescent eps = stochasticVanillaGradientAscent (-eps)+++--- Mean Squared Error --+++meanSquaredError+    :: ((c :#: m) -> [x] -> [Double]) -- ^ Error Function+    -> (c :#: m) -- ^ Current Point+    -> [x] -- ^ Sample points+    -> Double -- ^ Mean squared error+meanSquaredError err p xs =+    let rsdls = err p xs+     in (*0.5) . mean $ (**2) <$> rsdls+++--- Cauchy Sequences ---+++cauchyLimit :: Manifold m => Int -> Double -> [c :#: m] -> c :#: m+-- | Attempts to calculate the limit of a sequence. This finds the iterate with a+-- sufficiently small 'Euclidean' distance from the previous iterate, or returns+-- the nth iterate.+cauchyLimit n eps ps = last $ cauchySequence n eps ps++cauchySequence :: Manifold m => Int -> Double -> [c :#: m] -> [c :#: m]+cauchySequence n eps ps =+    let ps' = take n ps+        pps' = takeWhile taker . zip ps' $ tail ps'+     in head ps : map snd pps'+       where taker (p1,p2) = let p = alterChart Cartesian $ p1 <-> p2 in eps < sqrt (p <.> p)+++-- Least Squares --+++designMatrix :: Manifold m => [c :#: m] -> Function (Dual c) Cartesian :#: Tensor Euclidean m+-- | A glorified fromRows operation.+designMatrix rws = matrixTranspose $ coordinateTransform rws++leastSquares :: Manifold m => [c :#: m] -> [Double] -> Dual c :#: m+leastSquares xs ys =+    let mtx = designMatrix xs+        mtxt = matrixTranspose mtx+        prj = matrixInverse (mtxt <#> mtx) <#> mtxt+     in prj >.> euclideanPoint ys++leastSquares0 :: Manifold m => (Function c Cartesian :#: Tensor Euclidean m) -> [Double] -> c :#: m+leastSquares0 mtx ys =+    let mtxt = matrixTranspose mtx+        prj = matrixInverse (mtxt <#> mtx) <#> mtxt+     in prj >.> euclideanPoint ys++-- Newton --++newtonStep :: Manifold m+    => Double -- ^ Step size+    -> (Differentials :#: Tangent c m) -- ^ Derivatives+    -> (Function Partials Differentials :#: Tensor (Tangent c m) (Tangent c m)) -- ^ Hessian+    -> (c :#: m) -- ^ Step+newtonStep eps f' f'' = gradientStep (-eps) $ matrixInverse f'' >.> f'++newtonSequence :: Manifold m+    => Double -- ^ Step Size+    -> (c :#: m -> Differentials :#: Tangent c m) -- ^ Derivatives+    -> (c :#: m -> Function Partials Differentials :#: Tensor (Tangent c m) (Tangent c m)) -- ^ Hessian+    -> (c :#: m) -- ^ Initial point+    -> [c :#: m] -- ^ Newton sequence+newtonSequence eps f' f'' = iterate iterator+  where iterator p = newtonStep eps (f' p) (f'' p)+++-- Gauss Newton --++gaussNewtonStep :: Manifold m => Double -> [Double] -> [Differentials :#: Tangent c m] -> c :#: m+gaussNewtonStep eps rs grds = gradientStep (-eps) $ leastSquares0 (designMatrix grds) rs+++--- Graveyard ---+++{-+gaussNewtonPursuit :: Manifold m+    => Double -- ^ Damping Factor+    -> (c :#: m -> [x] -> [Double]) -- ^ Residual Function+    -> (c :#: m -> [x] -> [Differentials :#: Tangent c m]) -- ^ Residual Differential+    -> (c :#: m) -- ^ Initial guess+    -> StochasticPursuit x c m -- ^ Pursuit+gaussNewtonPursuit dmp residual residuald = accumulateFunction accumulator+  where accumulator xs cm =+            let cm' = gaussNewtonStep dmp (residual cm xs) (residuald cm xs)+             in (cm', cm')+             -}
+ Goal/Simulation/Physics/Configuration.hs view
@@ -0,0 +1,197 @@+-- | This module provides a general interface for working with mechanical systems.+module Goal.Simulation.Physics.Configuration+    ( -- * Configurations+    -- ** Types+      Generalized (Generalized)+    , PhaseSpace+    , GeneralizedVelocity+    , GeneralizedAcceleration+    , GeneralizedInertia+    -- ** Accessors+    , body+    , position+    , velocity+    , momentum+    -- * ForceFields+    -- ** Classes+    , ForceField (force)+    , Conservative (potentialEnergy)+    , vectorField+    , mechanicalEnergy+    -- ** Types+    , Gravity (Gravity)+    , earthGravity+    , Damping (Damping)+    -- * Util+    , periodic+    , revolutions+    , sliceVectorField+    , pairToPoint+    , pointToPair+    ) where++--- Imports ---+++-- Goal --++import Goal.Geometry++-- Qualified --++import qualified Data.Vector.Storable as C+++--- Configuration Space ---+++-- Charts --++-- | Generalized coordinates of a mechanical system.+data Generalized = Generalized++-- Manifolds --++-- | The 'Tangent' space of a mechanical system in 'Generalized' coordinates is the space of 'GeneralizedVelocity's.+type GeneralizedVelocity m = Tangent Generalized m++-- | The second order 'Tangent' space of a mechanical system is the space of 'GeneralizedAcceleration's.+type GeneralizedAcceleration m = Tangent Partials (GeneralizedVelocity m)++-- | The tangent 'Bundle' of a dynamical system is also known as the+-- 'PhaseSpace'. The "state" of a mechanical system is typically understood to+-- be an element of the 'PhaseSpace'.+type PhaseSpace m = Bundle Generalized m++-- | The 'Riemannian' metric on a mechanical system is known as the 'GeneralizedInertia'.+type GeneralizedInertia m = Tensor (GeneralizedVelocity m) (GeneralizedVelocity m)++-- Functions --++-- | Returns the underlying mechanical 'Manifold' of an element of the 'PhaseSpace'.+body :: Manifold m => Partials :#: PhaseSpace m -> m+body = removeBundle . manifold++-- | Returns the 'position' coordinates of a mechanical system.+position :: Manifold m => Partials :#: PhaseSpace m -> Generalized :#: m+position = projectTangent . bundleToTangent++-- | Returns the 'velocity' coordinates of a mechanical system.+velocity :: Manifold m => Partials :#: PhaseSpace m -> Partials :#: GeneralizedVelocity m+velocity = bundleToTangent++-- | Returns the 'momentum' coordinates of a mechanical system.+momentum :: Manifold m => Differentials :#: PhaseSpace m -> Differentials :#: GeneralizedVelocity m+momentum = bundleToTangent++kineticEnergy :: Riemannian Generalized m => Partials :#: GeneralizedVelocity m -> Double+kineticEnergy dq = 0.5 * (flat dq <.> dq)++-- Force fields --++-- | Gravitational force.+newtype Gravity = Gravity Double++earthGravity :: Gravity+earthGravity = Gravity 9.80665++newtype Damping = Damping Double++-- | A 'ForceField' describes how to attach a force vector (an element of the+-- dual space of generalized accelerations) to every point in the phase space.+-- Note that a 'ForceField' is not necessarily 'Conservative'.+class Manifold m => ForceField f m where+    force :: f -> Partials :#: PhaseSpace m -> Differentials :#: GeneralizedAcceleration m++-- | A 'Conservative' force depends only on 'position's and can be described as the gradient of a 'potentialEnergy'.+class ForceField f m => Conservative f m where+    potentialEnergy :: f -> Generalized :#: m -> Double++-- | The 'vectorField' function takes a 'ForceField' on a mechanical system and converts it into the appropriate element of the 'Tangent' space of the 'PhaseSpace'.+vectorField :: (Riemannian Generalized m, ForceField f m)+    => f+    -> Partials :#: PhaseSpace m+    -> Partials :#: Tangent Partials (PhaseSpace m)+vectorField f qdq = fromCoordinates (Tangent qdq) $ coordinates (velocity qdq) C.++ coordinates (sharp $ force f qdq)++mechanicalEnergy :: (Riemannian Generalized m, Conservative f m)+    => f+    -> Partials :#: PhaseSpace m+    -> (Double, Double)+mechanicalEnergy f qdq = (kineticEnergy $ velocity qdq, potentialEnergy f $ position qdq)+++--- Functions ---+++periodic :: Manifold m => [Bool] -> (Partials :#: PhaseSpace m) -> (Partials :#: PhaseSpace m)+periodic bls qdq =+    let q = position qdq+        q' = fromList (manifold q) $ clipper <$> zip bls (listCoordinates q)+        dq' = fromCoordinates (Tangent q') . coordinates $ velocity qdq+     in tangentToBundle dq'+       where clipper (True,x) = x - 2 * pi * fromIntegral (revolutions x)+             clipper (False,x) = x++revolutions :: Double -> Int+revolutions x+    | x >= pi = floor $ (x + pi) / (2*pi)+    | x <= -pi = ceiling $ (x - pi) / (2*pi)+    | otherwise = 0++sliceVectorField :: (Riemannian Generalized m, ForceField f m)+    => Double+    -> Int+    -> f+    -> Partials :#: PhaseSpace m+    -> (Double,Double)+    -> (Double,Double)+sliceVectorField scl i f qdq =+    pointToPair scl i . vectorField f . pairToPoint i qdq++pairToPoint :: Manifold m => Int -> Partials :#: PhaseSpace m -> (Double, Double) -> Partials :#: PhaseSpace m+pairToPoint n qdq (x,dx) =+    let (hxs,_:txs) = splitAt n . listCoordinates $ position qdq+        (hdxs,_:tdxs) = splitAt n . listCoordinates $ velocity qdq+     in fromList (manifold qdq) $ hxs ++ x:txs ++ hdxs ++ dx:tdxs++pointToPair :: Manifold m => Double -> Int -> (c :#: m) -> (Double, Double)+pointToPair scl n dqddq =+     (scl * coordinate n dqddq, scl * coordinate (n + div 2 (dimension $ manifold dqddq)) dqddq)++--- Instances ---+++instance (ForceField f m, ForceField g m) => ForceField (f,g) m where+    force (f,g) qdq = force f qdq  <+> force g qdq++instance Manifold m => ForceField (Partials :#: PhaseSpace m -> Differentials :#: GeneralizedAcceleration m) m where+    force f = f++-- Damping --++instance Manifold m => ForceField Damping m where+    force (Damping c) qdq =+        let dq = velocity qdq+         in fromCoordinates (Tangent dq) . coordinates $ alterCoordinates (negate . (*c)) dq+++--- Graveyard ---++{-+traversableToBundle :: (T.Traversable f, Manifold (f m), Manifold m) => f (d :#: Bundle c m) -> d :#: Bundle c (f m)+traversableToBundle qdqt =+    let pbt = Bundle $ removeBundle . manifold <$> qdqt+        (qs,qds) = unzip . F.toList $ cleanPoint <$> qdqt+     in fromCoordinates pbt $ C.concat qs C.++ C.concat qds+       where cleanPoint pbp = C.splitAt (dimension $ manifold pbp) $ coordinates pbp++bundleToTraversable :: (T.Traversable f, Manifold (f m), Manifold m) => d :#: Bundle c (f m) -> f (d :#: Bundle c m)+bundleToTraversable pbtp =+    let pbt = removeBundle $ manifold pbtp+     in snd . T.mapAccumL seperatePlanarLink (C.splitAt (dimension pbt) $ coordinates pbtp) $ pbt+      where seperatePlanarLink (qs,qds) pb =+                let (hqs,tqs) = C.splitAt (dimension pb) qs+                    (hqds,tqds) = C.splitAt (dimension pb) qds+                 in ((tqs,tqds), fromCoordinates (Bundle pb) $ hqs C.++ hqds)+-}
+ Goal/Simulation/Physics/Models/Pendulum.hs view
@@ -0,0 +1,41 @@+module Goal.Simulation.Physics.Models.Pendulum where+++--- Imports ---+++-- Goal --++import Goal.Geometry++import Goal.Simulation.Physics.Configuration+++--- Penduli ---+++data Pendulum = Pendulum Double Double deriving (Eq, Read, Show)+++--- Instances ---++instance Manifold Pendulum where+    dimension _ = 1++instance Riemannian Generalized Pendulum where+    metric q =+        let (Pendulum m l) = manifold q+         in fromList (Tensor (Tangent q) (Tangent q)) [m*l^2]++instance Conservative Gravity Pendulum where+    potentialEnergy (Gravity g) q =+        let (Pendulum m l) = manifold q+            [tht] = listCoordinates q+         in m * g * l * (1 - cos tht)++instance ForceField Gravity Pendulum where+    force (Gravity g) qdq =+        let q = position qdq+            (Pendulum m l) = manifold q+            [tht] = listCoordinates q+         in fromList (Tangent $ velocity qdq) [-m*g*l*sin tht]
+ Goal/Simulation/Plot.hs view
@@ -0,0 +1,109 @@+module Goal.Simulation.Plot+    ( module Goal.Simulation.Plot+    , mainGUI+    , initGUI )+    where+++-- Imports --+++-- Goal --++import Goal.Core hiding (on,set)+import Goal.Simulation.Mealy++-- Chart --++import Graphics.Rendering.Cairo as X hiding (lineTo,moveTo,x)+import Graphics.UI.Gtk++-- Unqualified --++import Data.IORef+import System.Clock+++--- Processes ---+++chainWindow :: Int -> Mealy x [x]+chainWindow n = accumulateMealy [] $ proc (x,xs) -> do+    let xs' = take n $ x:xs+    returnA -< (xs',xs')++trajectoryWindow :: Double -> Mealy (Double,x) [(Double,x)]+trajectoryWindow tivl = accumulateMealy [] $ proc ((t,x),xts) -> do+    let xts' = takeWhile (\(t',_) -> t - t' < tivl) $ (t,x):xts+    returnA -< (xts',xts')+++--- Animations ---+++--- IO ---+++data AnimationPost a = AP DrawingArea (Maybe Int) (IORef (Maybe (Renderable a))) (IORef TimeSpec)++changeFramerate :: Maybe Int -> AnimationPost a -> AnimationPost a+changeFramerate fps (AP da _ rnblrf tmrf) = AP da fps rnblrf tmrf++postRenderable :: AnimationPost a -> Renderable a -> IO ()+-- | Posts a renderable for animation. Note that this delays the calling thread until both+-- the image has been drawn and the frames per second interval has been passed. Also note+-- that this function is not in the business of frame skipping. (Although if I were to+-- use the difference information in both directions, it could be).+postRenderable (AP da Nothing rnblrf _) rnbl = do+    writeIORef rnblrf . Just $ rnbl+    postGUISync $ widgetQueueDraw da+postRenderable (AP da (Just fps) rnblrf tmrf) rnbl = do+    writeIORef rnblrf . Just $ rnbl+    postGUISync $ widgetQueueDraw da+    t0 <- readIORef tmrf+    t1 <- getTime Monotonic+    let diff = recip (fromIntegral fps) - fromIntegral (nsec t1 - nsec t0) / (10^9)+    threadDelay (round $ 10^6 * diff)+    writeIORef tmrf =<< getTime Monotonic++newAnimationPost ::  Double -> Maybe Int -> IO (AnimationPost a)+-- | Creates a new animation post, along with a window which can react to key+-- presses and within which the images will be drawn.+newAnimationPost art fps = do++    tmrf <- newIORef =<< getTime Monotonic+    rnblrf <- newIORef Nothing+    win <- windowNew+    afrm <- aspectFrameNew 0.5 0.5 . Just $ realToFrac art+    da <- drawingAreaNew++    void . (da `on` exposeEvent) . liftIO $ do++        rnbl <- readIORef rnblrf+        when (isJust rnbl) $ do+                void $ updateCanvas (fromJust rnbl) da+                return ()++        return True++{-+    (win `on` keyPressEvent) $ do++        ky <- eventKeyVal++        case keyToChar ky of+          Just cr ->++        liftIO $ print "BAM!"++        return True+        -}++    set afrm [ containerChild := da ]+    set win [ containerChild := afrm ]++    widgetShowAll win++    void $ onDestroy win mainQuit++    return $ AP da fps rnblrf tmrf
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2014, Sacha Sokoloski++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 Sacha Sokoloski 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
+ goal-simulation.cabal view
@@ -0,0 +1,120 @@+name: goal-simulation+version: 0.1+synopsis: Mealy based simulation tools+license: BSD3+license-file: LICENSE+author: Sacha Sokoloski+maintainer: sokolo@mis.mpg.de+category: Math+build-type: Simple+cabal-version: >=1.10++library+    exposed-modules:+        Goal.Simulation,+        Goal.Simulation.Mealy,+        Goal.Simulation.Chain,+        Goal.Simulation.Flow,+        Goal.Simulation.Filter,+        Goal.Simulation.Filter.Flow,+        Goal.Simulation.Optimization,+        Goal.Simulation.Plot,+        Goal.Simulation.Physics.Configuration,+        Goal.Simulation.Physics.Models.Pendulum+    default-extensions: TypeOperators, TypeFamilies, FlexibleInstances,+        FlexibleContexts, MultiParamTypeClasses, ScopedTypeVariables,+        RankNTypes, Arrows+    build-depends:+        base==4.*,+        goal-core==0.1,+        goal-geometry==0.1,+        goal-probability==0.1,+        machines==0.5.*,+        vector==0.11.*,+        hmatrix==0.17.*,+        cairo==0.13.*,+        gtk==0.14.*,+        clock==0.6.*+    default-language: Haskell2010+    ghc-options: -O2 -Wall -fno-warn-type-defaults -fno-warn-missing-signatures++executable rk4+    main-is: rk4.hs+    hs-source-dirs: scripts+    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults+        -fno-warn-missing-signatures -fno-warn-unused-do-bind+    build-depends: base, goal-core, goal-geometry, goal-simulation, vector+    default-language: Haskell2010++executable markov-chain+    main-is: markov-chain.hs+    hs-source-dirs: scripts+    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults+        -fno-warn-missing-signatures -fno-warn-unused-do-bind+    build-depends: base, goal-core, goal-geometry, goal-probability, vector,+        goal-simulation, hmatrix+    default-language: Haskell2010++executable ito-process+    main-is: ito-process.hs+    hs-source-dirs: scripts+    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults+        -fno-warn-missing-signatures -fno-warn-unused-do-bind+    build-depends: base, goal-core, goal-geometry, goal-probability, vector,+        goal-simulation, hmatrix+    default-language: Haskell2010++executable pendulum-vector-field+    main-is: vector-field.hs+    hs-source-dirs: scripts/pendulum+    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults+        -fno-warn-missing-signatures -fno-warn-unused-do-bind+    build-depends: base, goal-core, goal-geometry, goal-probability,+        goal-simulation, vector+    default-language: Haskell2010++executable pendulum-simulation+    main-is: simulation.hs+    hs-source-dirs: scripts/pendulum+    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults+        -fno-warn-missing-signatures -fno-warn-unused-do-bind+    build-depends: base, goal-core, goal-geometry, goal-probability,+        goal-simulation, vector+    default-language: Haskell2010++executable pendulum-filter-histogram+    main-is: filter-histogram.hs+    hs-source-dirs: scripts/pendulum+    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults+        -fno-warn-missing-signatures -fno-warn-unused-do-bind+    build-depends: base, goal-core, goal-geometry, goal-probability,+        goal-simulation, vector, directory+    default-language: Haskell2010++executable pendulum-filter-simulation+    main-is: filter-simulation.hs+    hs-source-dirs: scripts/pendulum+    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults+        -fno-warn-missing-signatures -fno-warn-unused-do-bind+    build-depends: base, goal-core, goal-geometry, goal-probability,+        goal-simulation, vector, directory+    default-language: Haskell2010++executable pendulum-filter-train+    main-is: filter-train.hs+    hs-source-dirs: scripts/pendulum+    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults+        -fno-warn-missing-signatures -fno-warn-unused-do-bind+    build-depends: base, goal-core, goal-geometry, goal-probability,+        goal-simulation, vector, directory+    default-language: Haskell2010++executable pendulum-filter-code+    main-is: filter-code.hs+    hs-source-dirs: scripts/pendulum+    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults+        -fno-warn-missing-signatures -fno-warn-unused-do-bind+    build-depends: base, goal-core, goal-geometry, goal-probability,+        goal-simulation, vector, directory, mtl+    default-language: Haskell2010+
+ scripts/ito-process.hs view
@@ -0,0 +1,52 @@+{-# LANGUAGE Arrows #-}++--- Imports ---+++-- Goal --++import Goal.Core++import Goal.Simulation+import Goal.Geometry+import Goal.Probability++import qualified Data.Vector.Storable as C+import qualified Numeric.LinearAlgebra.HMatrix as M++--- Main ---+++main :: IO ()+main = do+    let dt = 0.01+        fps = round $ recip dt+        tivl = 2+        t0 = 0+        ts = [t0,t0 + dt..]+        x0 = euclideanPoint [0]+    trj <- runWithSystemRandom $ itoProcess (\t _ -> C.fromList [cos t]) (\t _ -> M.fromLists [[1 + cos (2*t)]]) t0 x0++    let trajectoryToRenderable ln = toRenderable . execEC $ do++            layout_y_axis . laxis_generate .= scaledAxis def (-6,6)++            plot . liftEC $ do+                plot_lines_title .= "Ito Process"+                plot_lines_style .= solidLine 3 (opaque red)+                plot_lines_values .= [ln]++    let mly = proc t -> do+            x <- trj -< t+            pth <- trajectoryWindow tivl -< (t,coordinate 0 x)+            returnA -< trajectoryToRenderable pth++    --- Gtk ---++    initGUI++    apst <- newAnimationPost 2 (Just fps)++    forkIO $ streamM_ mly (postRenderable apst) ts++    mainGUI
+ scripts/markov-chain.hs view
@@ -0,0 +1,48 @@+{-# LANGUAGE Arrows #-}++--- Imports ---+++-- Goal --++import Goal.Core+import Goal.Geometry+import Goal.Simulation++import Goal.Probability++--- Main ---+++main :: IO ()+main = do+    let sz = 5+        nstps = 20+        ts = [1..] :: [Int]+        x0 = 1+        xs = [1..sz]+        smtx = fromList (markovTensor xs) . concat . replicate sz $ replicate sz (1/fromIntegral sz)++    chn <- runWithSystemRandom $ markovChain smtx x0++    let chainToRenderable ln = toRenderable . execEC $ do++           layout_y_axis . laxis_generate .= scaledIntAxis defaultIntAxis (0,sz + 1)++           plot . liftEC $ do+                   plot_lines_style .= solidLine 3 (opaque red)+                   plot_lines_title .= "Markov Chain"+                   plot_lines_values .= [ln]++    let mly = proc t -> do+            x <- chn -< ()+            stps <- chainWindow nstps -< (t,x)+            returnA -< chainToRenderable stps++    initGUI++    apst <- newAnimationPost 2 (Just 10)++    forkIO $ streamM_ mly (postRenderable apst) ts++    mainGUI
+ scripts/pendulum/filter-code.hs view
@@ -0,0 +1,87 @@++++-- Goal --++import Pendulum++import Goal.Core+import Goal.Geometry+import Goal.Probability++-- Qualified --++import qualified System.Directory as D+import qualified Control.Monad.State.Lazy as S+++--- Program ---++++-- Globals --++stp = 100++-- Functions --+++rateLayout clr mx rs = execEC $ do++    S.modify $ pixMapLayout nnrns nnrns++    layout_margin .= 1++    let clrs = [ clr `withOpacity` (r / mx) | r <- rs ]+        rcs = transpose $ reverse <$> breakEvery nnrns clrs++    layout_plots .= [pixMapPlot (0,0) rcs]++    plot . liftEC $ do++        let nnrns' = fromIntegral nnrns+        let vlns = [ [(k-0.5,-1.5),(k-0.5,nnrns' + 1)] | k <- [0..nnrns'] ]+            hlns = [ [(-1.5,k-0.5),(nnrns' + 1,k-0.5)] | k <- [0..nnrns'] ]++        plot_lines_style .= solidLine 1 (opaque grey)+        plot_lines_values .= (vlns ++ hlns)++-- Main --++main :: IO ()+main = do++    bl <- D.doesFileExist flnm+    c0s <- if bl+              then read <$> readFile flnm+              else error "Script requires a 'ppc-dynamics' file"++    let nnp = fromList nn c0s+        qdq0 = fromList (Bundle pndl) [2,0]++    xnzs <- runWithSystemRandom $ generatePath nstpssml nnp qdq0++    let xnz = xnzs !! stp+        xnz' = xnzs !! (stp+1)++    let (_,_,z) = splitTriple xnz+        (_,n',z') = splitTriple xnz'+        zs = listCoordinates z+        ns' = listCoordinates n'+        zs' = listCoordinates z'+        zs0' = zipWith (-) zs' ns'++        mx = maximum $ zs ++ zs'++        lytz = rateLayout red mx zs+        lytn' = rateLayout black mx ns'+        lytz0' = rateLayout blue mx zs0'+        lytz' = rateLayout red mx zs'+        rnbl = toRenderable . weights (1,1) $ tval lytz .|. tval lytz0' .|. tval lytn' .|. tval lytz'++    --void $ renderableToAspectWindow False 1200 300 rnbl+    putStrLn "Upper Bound:"+    print mx+    void $ renderableToFile (FileOptions (600,150) PDF) "filter-code.pdf" rnbl++
+ scripts/pendulum/filter-histogram.hs view
@@ -0,0 +1,112 @@+{-# LANGUAGE TypeOperators #-}++--- Imports ---+++-- Goal --++import Pendulum++import Goal.Core+import Goal.Geometry+import Goal.Probability+import Goal.Simulation++-- Qualified --++import qualified System.Directory as D+++--- Program ---+++-- Globals --++nbns = 20+mnbn = -2+mxbn = 6+npths = 1000+iaxprms = LinearAxisParams show 5 5+++-- Functions --++generatePaths+    :: Function Mixture Mixture :#: NeuralNetwork (Replicated Poisson) (Replicated Bernoulli) (Replicated Poisson)+    -> Rand (ST s) [[(Partials, Mixture, Mixture) :#: (PhaseSpace Pendulum, Replicated Poisson, Replicated Poisson)]]+generatePaths nnp = do+    let randomPath = randomState >>= generatePath nstpssml nnp+    replicateM npths randomPath++isLatentPathBounded+    :: [(Partials, Mixture, Mixture) :#: (PhaseSpace Pendulum, Replicated Poisson, Replicated Poisson)]+    -> Bool+isLatentPathBounded xnzs =+    let (xs,_,_) = unzip3 $ splitTriple <$> xnzs+     in and [ q > -pi && q < pi | [q,_] <- listCoordinates <$> xs ]++isLogLikelihoodPathFinite+    :: [(Partials, Mixture, Mixture) :#: (PhaseSpace Pendulum, Replicated Poisson, Replicated Poisson)]+    -> Bool+isLogLikelihoodPathFinite xnzs =+    let lzs = snd . unzip $ beliefNegativeLogLikelihoods trns <$> xnzs+     in and $ not . isInfinite <$> lzs++isLogLikelihoodPlotted+    :: (Partials, Mixture, Mixture) :#: (PhaseSpace Pendulum, Replicated Poisson, Replicated Poisson)+    -> Bool+isLogLikelihoodPlotted xnz =+    let lz = snd $ beliefNegativeLogLikelihoods trns xnz+     in lz > mnbn && lz < mxbn+++-- Main --++main :: IO ()+main = do++    bl <- D.doesFileExist flnm+    c0s <- if bl+              then read <$> readFile flnm+              else error "Script requires a 'ppc-dynamics' file"++    let nnp = fromList nn c0s++    xnzss <- runWithSystemRandom $ generatePaths nnp++    let bxnzss = filter isLatentPathBounded xnzss+        fbxnzss = filter isLogLikelihoodPathFinite bxnzss+        (pxnzs,pxnzs') = partition isLogLikelihoodPlotted $ concat fbxnzss+        (lns,lzs) = unzip $ beliefNegativeLogLikelihoods trns <$> pxnzs+        (lns',lzs') = unzip $ beliefNegativeLogLikelihoods trns <$> concat fbxnzss+        (_,ns,zs) = unzip3 $ splitTriple <$> pxnzs++    putStrLn "Percent in Bounds:"+    print $ 100 * genericLength bxnzss / genericLength xnzss+    putStrLn "Percent of Finite Rate:"+    print $ 100 * genericLength fbxnzss / genericLength bxnzss+    putStrLn "Percent within Histogram Bounds:"+    print $ 100 * genericLength pxnzs / genericLength (concat fbxnzss)+    putStrLn "Max out of Histogram Bounds:"+    print $ maximum . snd . unzip $ beliefNegativeLogLikelihoods trns <$> pxnzs'+    putStrLn "Average Observation Spike Count:"+    print . mean $ sum . listCoordinates <$> ns+    putStrLn "Average Belief Rate:"+    print . mean $ sum . listCoordinates <$> zs+    putStrLn "Average Information Gain:"+    print $ mean lns' - mean lzs'++    let plt =+            plot_bars_item_styles .~ [(FillStyleSolid $ opaque black,Nothing),(FillStyleSolid $ opaque red,Nothing)]+            $ histogramPlot nbns mnbn mxbn [lns,lzs] def++        lyt =+            layout_plots .~ [plotBars plt]+            -- $ layout_x_axis . laxis_title .~ "-Log-Likelihood"+            $ layout_y_axis . laxis_generate .~ autoScaledIntAxis iaxprms+            $ layout_x_axis . laxis_override .~ axisGridHide+            $ layout_y_axis . laxis_override .~ axisGridHide+            $ histogramLayout plt def++    --void $ renderableToAspectWindow False 1200 800 $ toRenderable (lyt :: Layout Double Int)+    void $ renderableToFile (FileOptions (600,200) PDF) "histogram.pdf" $ toRenderable (lyt :: Layout Double Int)
+ scripts/pendulum/filter-simulation.hs view
@@ -0,0 +1,109 @@++++-- Goal --++import Pendulum++import Goal.Core+import Goal.Geometry+import Goal.Probability+import Goal.Simulation++-- Qualified --++import qualified System.Directory as D+++--- Program ---+++-- Globals --++xaxprms = LinearAxisParams (show . round) 5 5+yaxprms = LinearAxisParams (show . round) 5 5++-- Functions --++coordinateLayout n rng (xs,ns,zs) = execEC $ do++    layout_x_axis . laxis_override .= axisGridHide+    layout_y_axis . laxis_generate .= scaledAxis yaxprms rng+    layout_y_axis . laxis_override .= axisGridHide+    layout_x_axis . laxis_generate .= autoScaledAxis xaxprms++    plot . liftEC $ do++        plot_lines_values .= [ zip ts $ coordinate n <$> xs ]+        plot_lines_style .= solidLine 2 (opaque black)++    plot . liftEC $ do++        plot_points_style .= filledCircles 2 (opaque black)+        plot_points_values .=+            zip ts (coordinate (2*n) . potentialMapping <$> (harmoniumTranspose trns >$> ns))++    plot . liftEC $ do++        plot_lines_style .= solidLine 2 (opaque red)+        plot_lines_values .=+            [zip ts $ coordinate (2*n) . potentialMapping <$> (harmoniumTranspose trns >$> zs)]+++-- Main --++main :: IO ()+main = do++    bl <- D.doesFileExist flnm+    c0s <- if bl+              then read <$> readFile flnm+              else error "Script requires a 'ppc-dynamics' file"++    let nnp = fromList nn c0s+        qdq0 = fromList (Bundle pndl) [1,0]++    xnzs <- runWithSystemRandom $ generatePath nstpssml nnp qdq0++    let (xs,ns,zs) = unzip3 $ splitTriple <$> xnzs+        (lns,lzs) = unzip $ beliefNegativeLogLikelihoods trns <$> xnzs++        gnlyt = execEC $ do++            layout_x_axis . laxis_override .= axisGridHide+            layout_y_axis . laxis_generate .= scaledAxis xaxprms (0,100)+            layout_y_axis . laxis_override .= axisGridHide+            layout_x_axis . laxis_generate .= autoScaledAxis xaxprms++            plot . liftEC $ do+                plot_points_values .= zip ts (sum . listCoordinates <$> ns)+                plot_points_style .= filledCircles 2 (opaque black)++            plot . liftEC $ do+                plot_lines_values .= [zip ts $ sum . listCoordinates <$> zs]+                plot_lines_style .= solidLine 2 (opaque red)++{-+            plot . liftEC $ do+                plot_lines_values .= [zip ts $ sum . listCoordinates <$> zipWith (<->) zs ns]+                plot_lines_style .= solidLine 2 (opaque blue)+                -}++        qlyt = coordinateLayout 0 (-2,2) (xs,ns,zs)+        dqlyt = coordinateLayout 1 (-4,4) (xs,ns,zs)++        rnbl = toRenderable $ StackedLayouts [StackedLayout qlyt, StackedLayout dqlyt, StackedLayout gnlyt] True++    putStrLn "Average Observation Likelihood:"+    print $ mean lns+    putStrLn "Average Belief Likelihood:"+    print $ mean lzs+    putStrLn "Average Observation Spike Count:"+    print . mean $ sum . listCoordinates <$> ns+    putStrLn "Average Belief Rate:"+    print . mean $ sum . listCoordinates <$> zs++    --void $ renderableToAspectWindow False 1200 1200 rnbl+    void $ renderableToFile (FileOptions (600,400) PDF) "simulation.pdf" rnbl++
+ scripts/pendulum/filter-train.hs view
@@ -0,0 +1,92 @@+{-# LANGUAGE BangPatterns,FlexibleContexts,DataKinds,TypeFamilies,TypeOperators,Arrows #-}++--- Imports ---+++-- Goal --++import Pendulum++import Goal.Core+import Goal.Geometry+import Goal.Probability+import Goal.Simulation++-- Qualified --++import qualified System.Directory as D+++--- Program ---+++-- Globals --++eps = 0.01+bnd = 2 :: Double+cdn = 5 :: Int+blkn = 10 :: Int+nstpstrn = 10 :: Int+nepchs = 100000 :: Int+nbns = 10++-- Functions --++trainerMWC nnp0 = do+    let randomPath nnp = randomState >>= generatePath nstpstrn nnp+    pathGenerator <- accumulateRandomFunction0 randomPath+    backpropagator <- boundedStochasticVanillaGradientDescent eps bnd (beliefBackpropagation trnss blkn cdn) nnp0+    return . accumulateMealy nnp0 $ proc ((),!nnp) -> do+        !nnp' <- backpropagator <<< pathGenerator -< nnp+        returnA -< (nnp',nnp')++-- Main --++main :: IO ()+main = do++    bl <- D.doesFileExist flnm+    c0s <- if bl+              then read <$> readFile flnm+              else runWithSystemRandom . replicateM (dimension nn) . generate . chart Standard $ fromList Normal [0,0.01]++    let nnp0 = fromList nn c0s++    trainer <- runWithSystemRandom $ trainerMWC nnp0++    --(lns,lzs,nnps) <- unzip3 <$> streamM trainer printerIO (replicate nepchs ())+    let nnp1 = last . take nepchs $ streamChain trainer+        (mp,mtx1,np,mtx2) = splitNeuralNetwork nnp1++    let wgtlyt = coordinateLogHistogram nbns "Network Weights" ["Second Layer Biases", "Second Layer", "First Layer Biases", "First Layer"]+            [coordinates mp, coordinates mtx1, coordinates np, coordinates mtx2]++    let qdq0 = fromList (Bundle pndl) [1.5,1.5]+        vflyt = execEC $ do++            layout_title .= "Vector Field"++            vectorFieldLayout++            plot . fmap plotVectorField . liftEC $ do++                vectorFieldPlot $ opaque black++                plot_vectors_title .= "Tru"+                plot_vectors_mapf .= sliceVectorField scl 0 f qdq0++            plot . fmap plotVectorField . liftEC $ do++                vectorFieldPlot $ opaque blue++                plot_vectors_mapf .= locationBeliefField trns0 dt scl 0 1 [0,0] nnp1+                plot_vectors_title .= "Est"+++    let rnbl = toRenderable $ StackedLayouts [StackedLayout vflyt, StackedLayout wgtlyt] False++    void $ renderableToAspectWindow False 400 800 rnbl++    writeFile flnm . show $ listCoordinates nnp1++
+ scripts/pendulum/simulation.hs view
@@ -0,0 +1,112 @@+{-# LANGUAGE TypeFamilies,FlexibleContexts,Arrows #-}++--- Imports ---+++-- Goal --++import Pendulum++import Goal.Core+import Goal.Geometry+import Goal.Probability+import Goal.Simulation+++--- Program ---+++-- Globals --++qdq0 = fromList (Bundle pndl) [1,0]++-- Functions --++pathToRenderable tqdqs =++    let phslyt = execEC $ do++            vectorFieldLayout++            layout_title .= "Phase Space"++            plot . liftEC $ do+                plot_lines_title .= "Phase"+                plot_lines_values .= [[(coordinate 0 qdq,coordinate 1 qdq) | (_,qdq) <- tqdqs]]+                plot_lines_style .= solidLine 3 (opaque black)++        enrlyt = execEC $ do++            layout_title .= "Energy"++            layout_x_axis . laxis_title .= "Time"++            layout_y_axis . laxis_generate .= scaledAxis def (0,20)+            layout_y_axis . laxis_title .= "Joules"++            let (tks,tus,tms) = unzip3 [ let (k,u) = mechanicalEnergy fg qdq in ((t,k),(t,u),(t,k+u)) | (t,qdq) <- tqdqs ]++            plot . liftEC $ do+                plot_lines_title .= "Kinetic"+                plot_lines_values .= [tks]+                plot_lines_style .= solidLine 3 (opaque red)++            plot . liftEC $ do+                plot_lines_title .= "Potential"+                plot_lines_values .= [tus]+                plot_lines_style .= solidLine 3 (opaque blue)++            plot . liftEC $ do+                plot_lines_title .= "Total"+                plot_lines_values .= [tms]+                plot_lines_style .= solidLine 3 (opaque purple)++        knmlyt = execEC $ do++            layout_title .= "Kinematics"++            let (_,qdq) = head tqdqs+                [tht] = listCoordinates $ position qdq+                (x,y) = (l * sin tht, -l * cos tht)++            layout_x_axis . laxis_generate .= scaledAxis def (-2,2)++            layout_y_axis . laxis_generate .= scaledAxis def (-2,2)++            plot . liftEC $ do+                plot_lines_values .= [[(0,0),(x,y)]]+                plot_lines_style .= solidLine 3 (opaque black)++            plot . liftEC $ do+                plot_points_values .= [(0,0)]+                plot_points_style .= hollowCircles 6 4 (opaque black)++            plot . liftEC $ do+                plot_points_values .= [(x,y)]+                plot_points_style .= filledCircles 6 (opaque black)+++     in toRenderable . weights (1,1) . wideAbove enrlyt $ tval knmlyt .|. tval phslyt++-- Main --++main :: IO ()+main = do++    lngvn <- runWithSystemRandom $ langevinFlow f sgma qdq0++    let tivl = 2.0+        fps = round $ recip dt++        mly = proc t -> do+            x <- lngvn -< t+            txs <- trajectoryWindow tivl -< (t,x)+            returnA -< pathToRenderable txs++    initGUI++    apst <- newAnimationPost 1.2 (Just fps)++    forkIO $ streamM_ mly (postRenderable apst) ts++    mainGUI
+ scripts/pendulum/vector-field.hs view
@@ -0,0 +1,46 @@+++--- Imports ---+++-- Goal --++import Pendulum++import Goal.Core+import Goal.Geometry+import Goal.Simulation+++--- Globals ---+++qdq0 = fromList (Bundle pndl) [0,0]+++--- Main ---+++main = do++    let rnbl = toRenderable . execEC $ do++            vectorFieldLayout++            layout_title .= "Pendulum Vector Field"++            plot . fmap plotVectorField . liftEC $ do++                vectorFieldPlot $ opaque red+                plot_vectors_mapf .= sliceVectorField 0.1 0 f qdq0+                plot_vectors_title .= "Damped"++            plot . fmap plotVectorField . liftEC $ do++                vectorFieldPlot $ opaque blue+                plot_vectors_mapf .= sliceVectorField 0.1 0 fg qdq0+                plot_vectors_title .= "Conservative"++    void $ renderableToAspectWindow False 800 800 rnbl+    --void $ renderableToFile (FileOptions (400,400) PDF) "vector-field.pdf" rnbl+
+ scripts/rk4.hs view
@@ -0,0 +1,123 @@+--- Imports ---+++-- Goal --++import Goal.Core+import Goal.Geometry++import Goal.Simulation++import qualified Data.Vector.Storable as C++--- Script ---+++main = do++    -- Generation --++    -- We can simulate sin either as non-autonomous or second order autonomous+    let sin' t _ = C.singleton $ cos t+        vsin' x = C.fromList [x C.! 1,-x C.! 0]+        exp' = id++        t0 = 0+        tf = 10+        dt1 = 2+        dt2 = 1+        dt3 = 0.1+        ts1 = [t0,t0+dt1..tf]+        ts2 = [t0,t0+dt2..tf]+        ts3 = [t0,t0+dt3..tf]++        sx0 = euclideanPoint [0]+        vx0 = euclideanPoint [0,1]+        ex0 = euclideanPoint [1]++        sinMealy = nonAutonomousODE sin' t0 sx0+        sinMealyEuler = nonAutonomousODEEuler sin' t0 sx0+        vsinMealy = autonomousODE vsin' vx0+        expMealy = autonomousODE exp' ex0+        expMealyEuler = autonomousODEEuler exp' ex0++        ssimulator ts mly = zip ts $ coordinate 0 <$> stream mly ts+        vsimulator ts mly = zip ts $ coordinate 0 <$> stream mly ts+        esimulator ts mly = zip ts $ coordinate 0 <$> stream mly ts++    -- Plots --++    -- Sin++    let sinrnbl = toRenderable . execEC $ do++            layout_title .= "Sin Wave (dt = {1,0.1})"++            plot . liftEC $ do+                plot_lines_style .= dashedLine 3 [10,5] (opaque black)+                plot_lines_title .= "True"+                plot_lines_values .= [zip ts3 $ sin <$> ts3]++            plot . liftEC $ do+                plot_lines_style .= solidLine 2 (opaque blue)+                plot_lines_title .= "RK4"+                plot_lines_values .=+                  [ ssimulator ts2 sinMealy, ssimulator ts3 sinMealy ]++            plot . liftEC $ do+                plot_lines_style .= solidLine 2 (opaque red)+                plot_lines_title .= "Euler"+                plot_lines_values .=+                  [ ssimulator ts2 sinMealyEuler, ssimulator ts3 sinMealyEuler ]++            plot . liftEC $ do+                plot_lines_style .= solidLine 2 (opaque purple)+                plot_lines_title .= "RK4 2nd Order"+                plot_lines_values .=+                  [ vsimulator ts2 vsinMealy , vsimulator ts3 vsinMealy ]++    -- Exponential++    let exprnbl = toRenderable . execEC $ do++            layout_title .= "Exponential Function (dt = {2,1,0.1})"++            plot . liftEC $ do+                plot_lines_style .= dashedLine 3 [10,5] (opaque black)+                plot_lines_title .= "True"+                plot_lines_values .= [zip ts3 $ exp <$> ts3]++            plot . liftEC $ do+                plot_lines_style .= solidLine 3 (opaque blue)+                plot_lines_title .= "RK4"+                plot_lines_values .= [ esimulator ts1 expMealy, esimulator ts2 expMealy, esimulator ts3 expMealy ]++            plot . liftEC $ do+                plot_lines_style .= solidLine 3 (opaque red)+                plot_lines_title .= "Euler"+                plot_lines_values .=+                    [ esimulator ts1 expMealyEuler, esimulator ts2 expMealyEuler, esimulator ts3 expMealyEuler ]++    -- IO++    renderableToAspectWindow False 800 600 . gridToRenderable . weights (1,1) . tallBeside sinrnbl $ tval exprnbl+++--- Extra Functions ---+++autonomousODEEuler f' p0 =+    accumulateFunction accumulator (0,p0)+      where accumulator t' (t,p) =+                let dt = t' - t+                    p' = p <+> fromCoordinates (manifold p) (stepEuler f' dt $ coordinates p)+                 in (p',(t',p'))++nonAutonomousODEEuler f' t0 p0 =+    accumulateFunction accumulator (t0,p0)+      where accumulator t' (t,p) =+                let dt = t' - t+                    p' = p <+> fromCoordinates (manifold p) (stepEuler' f' t dt $ coordinates p)+                 in (p',(t',p'))++