diff --git a/Goal/Simulation.hs b/Goal/Simulation.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Simulation.hs
@@ -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
diff --git a/Goal/Simulation/Chain.hs b/Goal/Simulation/Chain.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Simulation/Chain.hs
@@ -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')
diff --git a/Goal/Simulation/Filter.hs b/Goal/Simulation/Filter.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Simulation/Filter.hs
@@ -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')
+
diff --git a/Goal/Simulation/Filter/Flow.hs b/Goal/Simulation/Filter/Flow.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Simulation/Filter/Flow.hs
@@ -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)))
diff --git a/Goal/Simulation/Flow.hs b/Goal/Simulation/Flow.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Simulation/Flow.hs
@@ -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
diff --git a/Goal/Simulation/Mealy.hs b/Goal/Simulation/Mealy.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Simulation/Mealy.hs
@@ -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
diff --git a/Goal/Simulation/Optimization.hs b/Goal/Simulation/Optimization.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Simulation/Optimization.hs
@@ -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')
+             -}
diff --git a/Goal/Simulation/Physics/Configuration.hs b/Goal/Simulation/Physics/Configuration.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Simulation/Physics/Configuration.hs
@@ -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)
+-}
diff --git a/Goal/Simulation/Physics/Models/Pendulum.hs b/Goal/Simulation/Physics/Models/Pendulum.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Simulation/Physics/Models/Pendulum.hs
@@ -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]
diff --git a/Goal/Simulation/Plot.hs b/Goal/Simulation/Plot.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Simulation/Plot.hs
@@ -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
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -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.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/goal-simulation.cabal b/goal-simulation.cabal
new file mode 100644
--- /dev/null
+++ b/goal-simulation.cabal
@@ -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
+
diff --git a/scripts/ito-process.hs b/scripts/ito-process.hs
new file mode 100644
--- /dev/null
+++ b/scripts/ito-process.hs
@@ -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
diff --git a/scripts/markov-chain.hs b/scripts/markov-chain.hs
new file mode 100644
--- /dev/null
+++ b/scripts/markov-chain.hs
@@ -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
diff --git a/scripts/pendulum/filter-code.hs b/scripts/pendulum/filter-code.hs
new file mode 100644
--- /dev/null
+++ b/scripts/pendulum/filter-code.hs
@@ -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
+
+
diff --git a/scripts/pendulum/filter-histogram.hs b/scripts/pendulum/filter-histogram.hs
new file mode 100644
--- /dev/null
+++ b/scripts/pendulum/filter-histogram.hs
@@ -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)
diff --git a/scripts/pendulum/filter-simulation.hs b/scripts/pendulum/filter-simulation.hs
new file mode 100644
--- /dev/null
+++ b/scripts/pendulum/filter-simulation.hs
@@ -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
+
+
diff --git a/scripts/pendulum/filter-train.hs b/scripts/pendulum/filter-train.hs
new file mode 100644
--- /dev/null
+++ b/scripts/pendulum/filter-train.hs
@@ -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
+
+
diff --git a/scripts/pendulum/simulation.hs b/scripts/pendulum/simulation.hs
new file mode 100644
--- /dev/null
+++ b/scripts/pendulum/simulation.hs
@@ -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
diff --git a/scripts/pendulum/vector-field.hs b/scripts/pendulum/vector-field.hs
new file mode 100644
--- /dev/null
+++ b/scripts/pendulum/vector-field.hs
@@ -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
+
diff --git a/scripts/rk4.hs b/scripts/rk4.hs
new file mode 100644
--- /dev/null
+++ b/scripts/rk4.hs
@@ -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'))
+
+