diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -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 diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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. 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/app/Main.hs b/app/Main.hs new file mode 100644 --- /dev/null +++ b/app/Main.hs @@ -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 diff --git a/rhine-bayes.cabal b/rhine-bayes.cabal new file mode 100644 --- /dev/null +++ b/rhine-bayes.cabal @@ -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 diff --git a/src/Data/MonadicStreamFunction/Bayes.hs b/src/Data/MonadicStreamFunction/Bayes.hs new file mode 100644 --- /dev/null +++ b/src/Data/MonadicStreamFunction/Bayes.hs @@ -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 diff --git a/src/FRP/Rhine/Bayes.hs b/src/FRP/Rhine/Bayes.hs new file mode 100644 --- /dev/null +++ b/src/FRP/Rhine/Bayes.hs @@ -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