rhine-bayes (empty) → 0.8.1.1

raw patch · 8 files changed

+629/−0 lines, 8 filesdep +basedep +dunaidep +log-domainsetup-changed

Dependencies added: base, dunai, log-domain, mmorph, monad-bayes, rhine, rhine-bayes, rhine-gloss, time, transformers

Files

• ChangeLog.md +7/−0
• LICENSE +30/−0
• README.md +12/−0
• Setup.hs +2/−0
• app/Main.hs +351/−0
• rhine-bayes.cabal +89/−0
• src/Data/MonadicStreamFunction/Bayes.hs +53/−0
• src/FRP/Rhine/Bayes.hs +85/−0

@@ -0,0 +1,7 @@+# Revision history for rhine-gloss++## 0.8.1.1++* First version. Version numbers follow rhine.+* Introduces basic stochastic processes and Sequential Monte Carlo particle filter+* Thank you, Reuben Cohn-Gordon and Dominic Steinitz

@@ -0,0 +1,30 @@+Copyright (c) 2017, Manuel Bärenz++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 Manuel Bärenz 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.

@@ -0,0 +1,12 @@+# README++This package connects `rhine` to the [`monad-bayes`](hackage.haskell.org/package/monad-bayes) library for probabilistic programming and inference.+It provides:++* Some standard stochastic processes such as Brownian Motion and Levý processes+* A particle filter inference method called Sequential Monte Carlo++This allows you to do interactive probabilistic (i.e. involving randomness) programs,+and at the same time perform online inference, or realtime machine learning.+An example for this is given in `rhine-bayes/app/Main.hs`,+where inference is performed both on simulated values as well as external input given by the user.

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -0,0 +1,351 @@+{- | Interactive machine learning example.++In this example, you will find the following:++* A stochastic movement of a ball in heat bath,+  where the ball is kicked around in Brownian motion+* A noisy sensor observing the ball+* Particle filter inference trying to recover the actual (latent) position of the ball,+  as well as the temperature+* Visualization of the simulation and the inference result+* Two different architectures for the whole application:+  * A simple, noninteractive architecture where simulation, inference and visualization all run synchronously+  * A more scalable, modular, interactive architecture, where all these three systems run on separate clocks,+    and the user can interactively change the temperature of the heat bath+-}++-- base+import Control.Monad (void)+import Data.Maybe (fromMaybe)+import Data.Monoid (Product (Product, getProduct))+import GHC.Float (double2Float, float2Double)+import Text.Printf (printf)++-- transformers+import Control.Monad.Trans.Class++-- time+import Data.Time (addUTCTime, getCurrentTime)++-- mmorph+import Control.Monad.Morph++-- log-domain+import Numeric.Log hiding (sum)++-- monad-bayes+import Control.Monad.Bayes.Class hiding (posterior, prior)+import Control.Monad.Bayes.Population hiding (hoist)+import Control.Monad.Bayes.Sampler.Strict++-- rhine+import FRP.Rhine++-- rhine-gloss+import FRP.Rhine.Gloss.IO++-- rhine-bayes+import FRP.Rhine.Bayes+import FRP.Rhine.Gloss.Common++type Temperature = Double+type Pos = (Double, Double)+type Sensor = Pos++-- * Model++-- ** Prior++-- | Harmonic oscillator with white noise+prior1d ::+  (MonadDistribution m, Diff td ~ Double) =>+  -- | Starting position+  Double ->+  -- | Starting velocity+  Double ->+  BehaviourF m td Temperature Double+prior1d initialPosition initialVelocity = feedback 0 $ proc (temperature, position') -> do+  impulse <- arrM (normal 0) -< temperature+  let acceleration = (-3) * position' + impulse+  -- Integral over roughly the last 100 seconds, dying off exponentially, as to model a small friction term+  velocity <- arr (+ initialVelocity) <<< decayIntegral 10 -< acceleration+  position <- integralFrom initialPosition -< velocity+  returnA -< (position, position)++-- | 2D harmonic oscillator with noise+prior :: (MonadDistribution m, Diff td ~ Double) => BehaviourF m td Temperature Pos+prior = prior1d 10 0 &&& prior1d 0 10++-- ** Observation++-- | Internal utility because `gloss` operates on floats+double2FloatTuple :: (Double, Double) -> (Float, Float)+double2FloatTuple = double2Float *** double2Float++-- | An integral where the integrated value dies of exponentially+decayIntegral :: (VectorSpace v (Diff td), Monad m, Floating (Diff td)) => Diff td -> BehaviourF m td v v+decayIntegral timeConstant = (timeConstant *^) <$> average timeConstant++-- | The assumed standard deviation of the sensor noise+sensorNoiseTemperature :: Double+sensorNoiseTemperature = 1++-- | A generative model of the sensor noise+noise :: MonadDistribution m => Behaviour m td Pos+noise = whiteNoise sensorNoiseTemperature &&& whiteNoise sensorNoiseTemperature++-- | A generative model of the sensor position, given the noise+generativeModel :: (MonadDistribution m, Diff td ~ Double) => BehaviourF m td Pos Sensor+generativeModel = proc latent -> do+  noiseNow <- noise -< ()+  returnA -< latent ^+^ noiseNow++{- | This remodels the distribution defined by `noise` as a PDF,+   as to be used in the inference later.+-}+sensorLikelihood :: Pos -> Sensor -> Log Double+sensorLikelihood (posX, posY) (sensorX, sensorY) = normalPdf posX sensorNoiseTemperature sensorX * normalPdf posY sensorNoiseTemperature sensorY++-- ** User behaviour++-- | The initial value for the temperature, and also the initial guess for the temperature inference+initialTemperature :: Temperature+initialTemperature = 7++-- | We infer the temperature by randomly moving around with a Brownian motion (Wiener process).+temperatureProcess :: (MonadDistribution m, Diff td ~ Double) => BehaviourF m td () Temperature+temperatureProcess = proc () -> do+  temperatureFactor <- wienerLogDomain 20 -< ()+  returnA -< runLogDomain temperatureFactor * initialTemperature++-- | Auxiliary conversion function belonging to the log-domain library, see https://github.com/ekmett/log-domain/issues/38+runLogDomain :: Log Double -> Double+runLogDomain = exp . ln++-- * Filtering++{- | Generate a random position and sensor value, given a temperature.+   Used for simulating a situation upon which we will perform inference.+-}+genModelWithoutTemperature :: (MonadDistribution m, Diff td ~ Double) => BehaviourF m td Temperature (Sensor, Pos)+genModelWithoutTemperature = proc temperature -> do+  latent <- prior -< temperature+  sensor <- generativeModel -< latent+  returnA -< (sensor, latent)++{- | Given sensor data, sample a latent position and a temperature, and weight them according to the likelihood of the observed sensor position.+   Used to infer position and temperature.+-}+posteriorTemperatureProcess :: (MonadMeasure m, Diff td ~ Double) => BehaviourF m td Sensor (Temperature, Pos)+posteriorTemperatureProcess = proc sensor -> do+  temperature <- temperatureProcess -< ()+  latent <- prior -< temperature+  arrM score -< sensorLikelihood latent sensor+  returnA -< (temperature, latent)++-- | A collection of all displayable inference results+data Result = Result+  { temperature :: Temperature+  , measured :: Sensor+  , latent :: Pos+  , particles :: [((Temperature, Pos), Log Double)]+  }+  deriving (Show)++-- | The number of particles used in the filter. Change according to available computing power.+nParticles :: Int+nParticles = 100++-- * Visualization++{- | The monad in which our program will run.+   'SamplerIO' is for the probabilistic effects from @monad-bayes@,+   while 'GlossConcT' adds interactive effects from @gloss@.+-}+type App = GlossConcT SamplerIO++-- | Draw the results of the simulation and inference+visualisation :: Diff td ~ Double => BehaviourF App td Result ()+visualisation = proc Result{temperature, measured, latent, particles} -> do+  constMCl clearIO -< ()+  time <- sinceInitS -< ()+  arrMCl paintIO+    -<+      toThermometer $+        pictures+          [ translate 0 (-40) $ scale 0.2 0.2 $ color white $ pictures $ do+              (n, message) <-+                zip+                  [0 ..]+                  [ printf "Temperature: %.2f" temperature+                  , printf "Particles: %i" $ length particles+                  , printf "Time: %.1f" time+                  ]+              return $ translate 0 ((-150) * n) $ text message+          , color red $ rectangleUpperSolid thermometerWidth $ double2Float temperature * thermometerScale+          ]+  drawBall -< (measured, 0.3, red)+  drawBall -< (latent, 0.3, green)+  drawParticles -< particles++-- ** Parameters for the temperature display++thermometerPos :: (Float, Float)+thermometerPos = (-300, -300)++toThermometer :: Picture -> Picture+toThermometer = uncurry translate thermometerPos++thermometerScale :: Float+thermometerScale = 20++thermometerWidth :: Float+thermometerWidth = 20++-- ** Helpers for drawing elements of the visualization++drawBall :: BehaviourF App td (Pos, Double, Color) ()+drawBall = proc (position, width, theColor) -> do+  arrMCl paintIO -< scale 20 20 $ uncurry translate (double2FloatTuple position) $ color theColor $ circleSolid $ double2Float width++drawParticle :: BehaviourF App td ((Temperature, Pos), Log Double) ()+drawParticle = proc ((temperature, position), probability) -> do+  drawBall -< (position, 0.1, withAlpha (double2Float $ exp $ 0.2 * ln probability) white)+  arrMCl paintIO -< toThermometer $ translate 0 (double2Float temperature * thermometerScale) $ color (withAlpha (double2Float $ exp $ 0.2 * ln probability) white) $ rectangleSolid thermometerWidth 2++drawParticles :: BehaviourF App td [((Temperature, Pos), Log Double)] ()+drawParticles = proc particles -> do+  case particles of+    [] -> returnA -< ()+    p : ps -> do+      drawParticle -< p+      drawParticles -< ps++glossSettings :: GlossSettings+glossSettings =+  defaultSettings+    { display = InWindow "rhine-bayes" (1024, 960) (10, 10)+    }++-- * Integration++-- | There are different architectural choices for the whole application, these can be tested against each other+mains :: [(String, IO ())]+mains =+  [ ("single rate", mainSingleRate)+  , ("multi rate, temperature process", mainMultiRate)+  ]++main :: IO ()+main = do+  putStrLn $ ("Choose between: " ++) $ unwords $ zipWith (\n (title, _program) -> "\n" ++ show n ++ ": " ++ title) [1 ..] mains+  choice <- read <$> getLine+  map snd mains !! (choice - 1)++-- ** Single-rate : One simulation step = one inference step = one display step++{- | Given an actual temperature, simulate a latent position and measured sensor position,+   and based on the sensor data infer the latent position and the temperature.+-}+filtered :: Diff td ~ Double => BehaviourF App td Temperature Result+filtered = proc temperature -> do+  (measured, latent) <- genModelWithoutTemperature -< temperature+  particles <- runPopulationCl nParticles resampleSystematic posteriorTemperatureProcess -< measured+  returnA+    -<+      Result+        { temperature+        , measured+        , latent+        , particles+        }++-- | Run simulation, inference, and visualization synchronously+mainClSF :: Diff td ~ Double => BehaviourF App td () ()+mainClSF = proc () -> do+  output <- filtered -< initialTemperature+  visualisation -< output++-- | Rescale to the 'Double' time domain+type GlossClock = RescaledClock GlossSimClockIO Double++glossClock :: GlossClock+glossClock =+  RescaledClock+    { unscaledClock = GlossSimClockIO+    , rescale = float2Double+    }++mainSingleRate =+  void $+    sampleIO $+      launchGlossThread glossSettings $+        reactimateCl glossClock mainClSF++-- ** Multi-rate: Simulation, inference, display at different rates++-- | Rescale the gloss clocks so they will be compatible with real 'UTCTime' (needed for compatibility with 'Millisecond')+type GlossClockUTC cl = RescaledClockS (GlossConcT IO) cl UTCTime (Tag cl)++glossClockUTC :: Real (Time cl) => cl -> GlossClockUTC cl+glossClockUTC cl =+  RescaledClockS+    { unscaledClockS = cl+    , rescaleS = const $ do+        now <- liftIO getCurrentTime+        return (arr $ \(timePassed, event) -> (addUTCTime (realToFrac timePassed) now, event), now)+    }++{- | The part of the program which simulates latent position and sensor,+   running 100 times a second.+-}+modelRhine :: Rhine (GlossConcT IO) (LiftClock IO GlossConcT (Millisecond 100)) Temperature (Temperature, (Sensor, Pos))+modelRhine = hoistClSF sampleIOGloss (clId &&& genModelWithoutTemperature) @@ liftClock waitClock++-- | The user can change the temperature by pressing the up and down arrow keys.+userTemperature :: ClSF (GlossConcT IO) (GlossClockUTC GlossEventClockIO) () Temperature+userTemperature = tagS >>> arr (selector >>> fmap Product) >>> mappendS >>> arr (fmap getProduct >>> fromMaybe 1 >>> (* initialTemperature))+ where+  selector (EventKey (SpecialKey KeyUp) Down _ _) = Just 1.2+  selector (EventKey (SpecialKey KeyDown) Down _ _) = Just (1 / 1.2)+  selector _ = Nothing++{- | This part performs the inference (and passes along temperature, sensor and position simulations).+   It runs as fast as possible, so this will potentially drain the CPU.+-}+inference :: Rhine (GlossConcT IO) (LiftClock IO GlossConcT Busy) (Temperature, (Sensor, Pos)) Result+inference = hoistClSF sampleIOGloss inferenceBehaviour @@ liftClock Busy+ where+  inferenceBehaviour :: (MonadDistribution m, Diff td ~ Double, MonadIO m) => BehaviourF m td (Temperature, (Sensor, Pos)) Result+  inferenceBehaviour = proc (temperature, (measured, latent)) -> do+    particles <- runPopulationCl nParticles resampleSystematic posteriorTemperatureProcess -< measured+    returnA -< Result{temperature, measured, latent, particles}++-- | Visualize the current 'Result' at a rate controlled by the @gloss@ backend, usually 30 FPS.+visualisationRhine :: Rhine (GlossConcT IO) (GlossClockUTC GlossSimClockIO) Result ()+visualisationRhine = hoistClSF sampleIOGloss visualisation @@ glossClockUTC GlossSimClockIO++-- | Compose all four asynchronous components to a single 'Rhine'.+mainRhineMultiRate =+  userTemperature+    @@ glossClockUTC GlossEventClockIO+      >-- keepLast initialTemperature -@- glossConcurrently -->+        modelRhine+        >-- keepLast (initialTemperature, (zeroVector, zeroVector)) -@- glossConcurrently -->+          inference+            >-- keepLast Result{temperature = initialTemperature, measured = zeroVector, latent = zeroVector, particles = []} -@- glossConcurrently -->+              visualisationRhine++mainMultiRate :: IO ()+mainMultiRate =+  void $+    launchGlossThread glossSettings $+      flow mainRhineMultiRate++-- * Utilities++instance MonadDistribution m => MonadDistribution (GlossConcT m) where+  random = lift random++sampleIOGloss :: App a -> GlossConcT IO a+sampleIOGloss = hoist sampleIO

@@ -0,0 +1,89 @@+name:                rhine-bayes+version:             0.8.1.1+synopsis:            monad-bayes backend for Rhine+description:+  This package provides a backend to the `monad-bayes` library,+  enabling you to write stochastic processes as signal functions,+  and performing online machine learning on them.+license:             BSD3+license-file:        LICENSE+author:              Manuel Bärenz+maintainer:          programming@manuelbaerenz.de+-- copyright:+category:            FRP+build-type:          Simple+extra-source-files:  ChangeLog.md+extra-doc-files:     README.md+cabal-version:       1.18++source-repository head+  type:     git+  location: git@github.com:turion/rhine.git++source-repository this+  type:     git+  location: git@github.com:turion/rhine.git+  tag:      v0.8.1.1++library+  exposed-modules:+    FRP.Rhine.Bayes+  other-modules:+    Data.MonadicStreamFunction.Bayes+  build-depends:       base         >= 4.11 && < 4.18+                     , transformers >= 0.5+                     , rhine        == 0.8.1.1+                     , dunai        >= 0.8+                     , log-domain   >= 0.12+                     , monad-bayes  >= 1.1.0+  hs-source-dirs:      src+  default-language:    Haskell2010+  default-extensions:+    Arrows+    DataKinds+    DeriveFunctor+    FlexibleContexts+    FlexibleInstances+    GeneralizedNewtypeDeriving+    MultiParamTypeClasses+    RankNTypes+    ScopedTypeVariables+    TupleSections+    TypeFamilies++  ghc-options:         -W+  if flag(dev)+    ghc-options: -Werror++executable rhine-bayes-gloss+  main-is:             Main.hs+  hs-source-dirs:      app+  ghc-options:         -threaded+  build-depends:       base         >= 4.11 && < 4.18+                     , rhine+                     , rhine-bayes+                     , rhine-gloss+                     , monad-bayes+                     , transformers+                     , log-domain+                     , mmorph+                     , time+  default-language:    Haskell2010+  default-extensions:+    Arrows+    DataKinds+    FlexibleContexts+    NamedFieldPuns+    RankNTypes+    TupleSections+    TypeApplications+    TypeFamilies++  ghc-options:         -W -threaded -rtsopts -with-rtsopts=-N+  if flag(dev)+    ghc-options: -Werror++flag dev+  description: Enable warnings as errors. Active on ci.+  default: False+  manual: True

@@ -0,0 +1,53 @@+module Data.MonadicStreamFunction.Bayes where++-- base+import Control.Arrow+import Data.Functor (($>))+import Data.Tuple (swap)++-- transformers++-- log-domain+import Numeric.Log hiding (sum)++-- monad-bayes+import Control.Monad.Bayes.Population++-- dunai+import Data.MonadicStreamFunction+import Data.MonadicStreamFunction.InternalCore (MSF (..))++-- | Run the Sequential Monte Carlo algorithm continuously on an 'MSF'+runPopulationS ::+  forall m a b.+  Monad m =>+  -- | Number of particles+  Int ->+  -- | Resampler+  (forall x. Population m x -> Population m x) ->+  MSF (Population m) a b ->+  -- FIXME Why not MSF m a (Population b)+  MSF m a [(b, Log Double)]+runPopulationS nParticles resampler = runPopulationsS resampler . (spawn nParticles $>)++-- | Run the Sequential Monte Carlo algorithm continuously on a 'Population' of 'MSF's+runPopulationsS ::+  Monad m =>+  -- | Resampler+  (forall x. Population m x -> Population m x) ->+  Population m (MSF (Population m) a b) ->+  MSF m a [(b, Log Double)]+runPopulationsS resampler = go+ where+  go msfs = MSF $ \a -> do+    -- TODO This is quite different than the dunai version now. Maybe it's right nevertheless.+    -- FIXME This normalizes, which introduces bias, whatever that means+    bAndMSFs <- runPopulation $ normalize $ resampler $ flip unMSF a =<< msfs+    return $+      second (go . fromWeightedList . return) $+        unzip $+          (swap . fmap fst &&& swap . fmap snd) . swap <$> bAndMSFs++-- FIXME see PR re-adding this to monad-bayes+normalize :: Monad m => Population m a -> Population m a+normalize = fromWeightedList . fmap (\particles -> second (/ (sum $ snd <$> particles)) <$> particles) . runPopulation

@@ -0,0 +1,85 @@+module FRP.Rhine.Bayes where++-- log-domain+import Numeric.Log hiding (sum)++-- monad-bayes+import Control.Monad.Bayes.Class+import Control.Monad.Bayes.Population++-- dunai+import qualified Control.Monad.Trans.MSF.Reader as DunaiReader++-- dunai-bayes+import qualified Data.MonadicStreamFunction.Bayes as DunaiBayes++-- rhine+import FRP.Rhine++-- * Inference methods++-- | Run the Sequential Monte Carlo algorithm continuously on a 'ClSF'.+runPopulationCl :: forall m cl a b . Monad m =>+  -- | Number of particles+  Int ->+  -- | Resampler (see 'Control.Monad.Bayes.Population' for some standard choices)+  (forall x . Population m x -> Population m x)+  -- | A signal function modelling the stochastic process on which to perform inference.+  --   @a@ represents observations upon which the model should condition, using e.g. 'score'.+  --   It can also additionally contain hyperparameters.+  --   @b@ is the type of estimated current state.+  -> ClSF (Population m) cl a b+  -> ClSF m cl a [(b, Log Double)]+runPopulationCl nParticles resampler = DunaiReader.readerS . DunaiBayes.runPopulationS nParticles resampler . DunaiReader.runReaderS++-- * Short standard library of stochastic processes++-- | White noise, that is, an independent normal distribution at every time step.+whiteNoise :: MonadDistribution m => Double -> Behaviour m td Double+whiteNoise sigma = constMCl $ normal 0 sigma++-- | Construct a Lévy process from the increment between time steps.+levy ::+  (MonadDistribution m, VectorSpace v (Diff td)) =>+  -- | The increment function at every time step. The argument is the difference between times.+  (Diff td -> m v) ->+  Behaviour m td v+levy incrementor = sinceLastS >>> arrMCl incrementor >>> sumS++-- | The Wiener process, also known as Brownian motion.+wiener, brownianMotion ::+  (MonadDistribution m, Diff td ~ Double) =>+  -- | Time scale of variance.+  Diff td ->+  Behaviour m td Double+wiener timescale = levy $ \diffTime -> normal 0 $ sqrt $ diffTime / timescale++brownianMotion = wiener++-- | The Wiener process, also known as Brownian motion, with varying variance parameter.+wienerVarying, brownianMotionVarying ::+  (MonadDistribution m, Diff td ~ Double) =>+  BehaviourF m td (Diff td) Double+wienerVarying = proc timeScale -> do+  diffTime <- sinceLastS -< ()+  let stdDev = sqrt $ diffTime / timeScale+  increment <- if stdDev > 0+    then arrM $ normal 0 -< stdDev+    else returnA -< 0+  sumS -< increment++brownianMotionVarying = wienerVarying++-- | The 'wiener' process transformed to the Log domain, also called the geometric Wiener process.+wienerLogDomain ::+  (MonadDistribution m, Diff td ~ Double) =>+  -- | Time scale of variance+  Diff td ->+  Behaviour m td (Log Double)+wienerLogDomain timescale = wiener timescale >>> arr Exp++-- | See 'wienerLogDomain' and 'wienerVarying'.+wienerVaryingLogDomain ::+  (MonadDistribution m, Diff td ~ Double) =>+  BehaviourF m td (Diff td) (Log Double)+wienerVaryingLogDomain = wienerVarying >>> arr Exp