flat-mcmc 0.1.0.0 → 1.0.0

raw patch · 7 files changed

+322/−200 lines, 7 filesdep +flat-mcmcdep +mcmc-typesdep +mwc-probabilitydep −mtldep −mwc-randomdep ~basedep ~monad-pardep ~monad-par-extrassetup-changedPVP ok

version bump matches the API change (PVP)

Dependencies added: flat-mcmc, mcmc-types, mwc-probability, pipes, transformers

Dependencies removed: mtl, mwc-random

Dependency ranges changed: base, monad-par, monad-par-extras, primitive, vector

API changes (from Hackage documentation)

Files

• LICENSE +16/−26
• Numeric/MCMC/Flat.hs +0/−157
• Setup.hs +5/−0
• flat-mcmc.cabal +74/−17
• lib/Numeric/MCMC/Flat.hs +183/−0
• test/BNN.hs +22/−0
• test/Rosenbrock.hs +22/−0

@@ -1,29 +1,19 @@-Copyright (C) 2012 Jared Tobin--Redistribution and use in source and binary forms, with or without-modification, are permitted provided that the following conditions are-met:--   1. Redistributions of source code must retain the above copyright-      notice, this list of conditions, and the following disclaimer.+Copyright (c) 2012-2016 Jared Tobin -   2. 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.+Permission is hereby granted, free of charge, to any person obtaining a copy+of this software and associated documentation files (the "Software"), to deal+in the Software without restriction, including without limitation the rights+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell+copies of the Software, and to permit persons to whom the Software is+furnished to do so, subject to the following conditions: -   3. The name of the author may not be used to endorse or promote-      products derived from this software without specific prior-      written permission.+The above copyright notice and this permission notice shall be included in+all copies or substantial portions of the Software. -THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``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 AUTHOR 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.+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN+THE SOFTWARE.

@@ -1,157 +0,0 @@-{-# OPTIONS_GHC -Wall #-}-{-# LANGUAGE BangPatterns #-}--module Numeric.MCMC.Flat (-            MarkovChain(..), Options(..), Ensemble-          , runChain, readInits-          ) where--import Control.Arrow-import Control.Monad-import Control.Monad.Reader-import Control.Monad.Primitive-import System.Random.MWC-import qualified Data.Vector           as V-import qualified Data.Vector.Unboxed   as U-import Control.Monad.Par                    (NFData)-import Control.Monad.Par.Scheds.Direct-import Control.Monad.Par.Combinator-import GHC.Float-import System.IO---- | Parallel map with a specified granularity.-parMapChunk :: NFData b => Int -> (a -> b) -> [a] -> Par [b]-parMapChunk n f xs = do-    xss <- parMap (map f) (chunk n xs)-    return (concat xss)-  where chunk _ [] = []-        chunk m ys = let (as, bs) = splitAt m ys-                     in  as : chunk m bs---- | State of the Markov chain.  Current ensemble position is held in 'theta',---   while 'accepts' counts the number of proposals accepted.-data MarkovChain = MarkovChain { ensemble :: Ensemble      -                               , accepts  :: {-# UNPACK #-} !Int }---- | Display the current state.  This will be very slow and should be replaced.-instance Show MarkovChain where-    show config = filter (`notElem` "[]") $ unlines $ map (show . map double2Float) (V.toList (ensemble config))---- | Options for the chain.  The target (expected to be a log density), as---   well as the size of the ensemble.  The size should be an even number.  Also---   holds the specified parallel granularity as 'csize'.-data Options = Options { _target :: [Double] -> Double    -                       , _size   :: {-# UNPACK #-} !Int -                       , _csize  :: {-# UNPACK #-} !Int } ---- | An ensemble of particles.-type Ensemble = V.Vector [Double]---- | A result with this type has a view of the chain's options.-type ViewsOptions = ReaderT Options---- | Generate a random value from a distribution having the property that ---   g(1/z) = zg(z).-symmetricVariate :: PrimMonad m => Gen (PrimState m) -> m Double-symmetricVariate g = do-    z <- uniformR (0 :: Double, 1 :: Double) g-    return $! 0.5*(z + 1)^(2 :: Int)-{-# INLINE symmetricVariate #-}---- | The result of a single-particle Metropolis accept/reject step.  This ---   compares a particle's state to a perturbation made by an affine ---   transformation based on a complementary particle.  Non-monadic to ---   more easily be used in the Par monad.-metropolisResult :: [Double] -> [Double] -- Target and alternate particles-                 -> Double -> Double     -- z ~ g(z) and zc ~ rand-                 -> ([Double] -> Double) -- Target function-                 -> ([Double], Int)      -- Result and accept counter-metropolisResult w0 w1 z zc target = -    let val      = target proposal - target w0 + (fromIntegral (length w0) - 1) * log z-        proposal = zipWith (+) (map (*z) w0) (map (*(1-z)) w1) -    in  if zc <= min 1 (exp val) then (proposal, 1) else (w0, 0)-{-# INLINE metropolisResult #-}---- | Execute Metropolis steps on the particles of a sub-ensemble by---   perturbing them with affine transformations based on particles---   in a complementary ensemble, in parallel.-executeMoves :: (Functor m, PrimMonad m)-             => Ensemble                            -- Target sub-ensemble-             -> Ensemble                            -- Complementary sub-ensemble-             -> Int                                 -- Size of the sub-ensembles-             -> Gen (PrimState m)                   -- MWC PRNG-             -> ViewsOptions m (Ensemble, Int)      -- Updated ensemble and # of accepts-executeMoves e0 e1 n g = do-    Options t _ csize <- ask--    zs  <- replicateM n (lift $ symmetricVariate g)-    zcs <- replicateM n (lift $ uniformR (0 :: Double, 1 :: Double) g)-    js  <- fmap U.fromList (replicateM n (lift $ uniformR (1:: Int, n) g))--    let w0 k    = e0 `V.unsafeIndex` (k - 1) -        w1 k ks = e1 `V.unsafeIndex` ((ks `U.unsafeIndex` (k - 1)) - 1) --        result  = runPar $ parMapChunk csize -                     (\(k, z, zc) -> metropolisResult (w0 k) (w1 k js) z zc t) -                         (zip3 [1..n] zs zcs)-        (newstate, nacc) = (V.fromList . map fst &&& sum . map snd) result--    return (newstate, nacc)-{-# INLINE executeMoves #-}---- | Perform a Metropolis accept/reject step on the ensemble by---   perturbing each element and accepting/rejecting the perturbation in---   parallel.-metropolisStep :: (Functor m, PrimMonad m)-               => MarkovChain                   -- State of the Markov chain-               -> Gen (PrimState m)             -- MWC PRNG-               -> ViewsOptions m MarkovChain    -- Updated sub-ensemble-metropolisStep state g = do-    Options _ n _ <- ask-    let n0        = truncate (fromIntegral n / (2 :: Double)) :: Int-        (e, nacc) = (ensemble &&& accepts) state-        (e0, e1)  = (V.slice (0 :: Int) n0 &&& V.slice n0 n0) e- -    -- Update each sub-ensemble -    result0 <- executeMoves e0 e1            n0 g-    result1 <- executeMoves e1 (fst result0) n0 g--    return $! -      MarkovChain (V.concat $ map fst [result0, result1]) -                  (nacc + snd result0 + snd result1)-{-# INLINE metropolisStep #-}---- | Diffuse through states.-runChain :: Options         -- Options of the Markov chain-         -> Int             -- Number of epochs to iterate the chain-         -> Int             -- Print every nth iteration.-         -> MarkovChain     -- Initial state of the Markov chain-         -> Gen RealWorld   -- MWC PRNG-         -> IO MarkovChain  -- End state of the Markov chain, wrapped in IO-runChain opts nepochs thinEvery initConfig g -    | l == 0 -        = error "runChain: ensemble must contain at least one particle"-    | l < (length . V.head) (ensemble initConfig)-        = do hPutStrLn stderr $ "runChain: ensemble should be twice as large as "-                             ++ "the target's dimension.  Continuing anyway."-             go opts nepochs thinEvery initConfig g-    | otherwise = go opts nepochs thinEvery initConfig g-  where -    l = V.length (ensemble initConfig)-    go o n t !c g0 | n == 0 = return c-                   | n `rem` t /= 0 = do-                       r <- runReaderT (metropolisStep c g0) o-                       go o (n - 1) t r g0-                   | otherwise = do-                       r <- runReaderT (metropolisStep c g0) o-                       print r-                       go o (n - 1) t r g0-{-# INLINE runChain #-}---- | A convenience function to read and parse ensemble inits from disk.  ---   Assumes a text file with one particle per line, where each particle---   element is separated by whitespace.-readInits :: FilePath -> IO Ensemble-readInits p = fmap (V.fromList . map (map read . words) . lines) (readFile p)-{-# INLINE readInits #-}-

@@ -1,2 +1,7 @@+module Main (main) where+ import Distribution.Simple++main :: IO () main = defaultMain+

@@ -1,31 +1,88 @@--- Initial flat_mcmc.cabal generated by cabal init.  For further --- documentation, see http://haskell.org/cabal/users-guide/- name:                flat-mcmc-version:             0.1.0.0+version:             1.0.0 synopsis:            Painless general-purpose sampling.--- description:          homepage:            http://jtobin.github.com/flat-mcmc--- license:             -license:             BSD3+license:             MIT license-file:        LICENSE author:              Jared Tobin maintainer:          jared@jtobin.ca--- copyright:           -category:            Numeric, Machine Learning, Statistics+category:            Math build-type:          Simple-cabal-version:       >=1.8-Description:--    Painless general-purpose sampling.+cabal-version:       >=1.10+description:+  *flat-mcmc* is a Haskell library for painless, efficient, general-purpose+  sampling from continuous distributions.+  .+  *flat-mcmc* uses an ensemble sampler that is invariant to affine+  transformations of space.  It wanders a target probability distribution's+  parameter space as if it had been "flattened" or "unstretched" in some sense,+  allowing many particles to explore it locally and in parallel.+  .+  In general this sampler is useful when you want decent performance without+  dealing with any tuning parameters or local proposal distributions.+  .+  *flat-mcmc* exports an 'mcmc' function that prints a trace to stdout, as well+  as a 'flat' transition operator that can be used more generally.+  .+  > import Numeric.MCMC.Flat+  > import Data.Vector (Vector, toList, fromList)+  >+  > rosenbrock :: Vector Double -> Double+  > rosenbrock xs = negate (5  *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2) where+  >   [x0, x1] = toList xs+  >+  > ensemble :: Ensemble+  > ensemble = fromList [+  >     fromList [negate 1.0, negate 1.0]+  >   , fromList [negate 1.0, 1.0]+  >   , fromList [1.0, negate 1.0]+  >   , fromList [1.0, 1.0]+  >   ]+  >+  > main :: IO ()+  > main = withSystemRandom . asGenIO $ mcmc 12500 ensemble rosenbrock  Source-repository head   Type:     git   Location: http://github.com/jtobin/flat-mcmc.git  library-  exposed-modules:     Numeric.MCMC.Flat-  -- other-modules:       -  build-depends:       base ==4.*, mtl ==2.1.*, primitive ==0.4.*, mwc-random ==0.12.*, vector ==0.9.*, monad-par ==0.3.*, monad-par-extras ==0.3.*-  ghc-options:         -Wall+  default-language: Haskell2010+  ghc-options:      -Wall+  hs-source-dirs:   lib+  exposed-modules:  Numeric.MCMC.Flat+  build-depends:+      base             <  5+    , mcmc-types       >= 1.0.1 && < 2+    , monad-par+    , monad-par-extras+    , mwc-probability  >= 1.0.1 && < 2+    , pipes            >= 4 && < 5+    , primitive+    , transformers+    , vector++Test-suite rosenbrock+  type:                exitcode-stdio-1.0+  hs-source-dirs:      test+  main-is:             Rosenbrock.hs+  default-language:    Haskell2010+  ghc-options:+    -rtsopts -threaded+  build-depends:+      base              < 5+    , flat-mcmc+    , vector++Test-suite bnn+  type:                exitcode-stdio-1.0+  hs-source-dirs:      test+  main-is:             BNN.hs+  default-language:    Haskell2010+  ghc-options:+    -rtsopts -threaded+  build-depends:+      base              < 5+    , flat-mcmc+    , vector 

@@ -0,0 +1,183 @@+{-# OPTIONS_GHC -fno-warn-type-defaults #-}+{-# LANGUAGE RecordWildCards #-}++-- |+-- Module: Numeric.MCMC.Flat+-- Copyright: (c) 2016 Jared Tobin+-- License: MIT+--+-- Maintainer: Jared Tobin <jared@jtobin.ca>+-- Stability: unstable+-- Portability: ghc+--+-- This is the 'affine invariant ensemble' or AIEMCMC algorithm described in+-- Goodman and Weare, 2010.  It is a reasonably efficient and hassle-free+-- sampler, requiring no mucking with tuning parameters or local proposal+-- distributions.+--+-- The 'mcmc' function streams a trace to stdout to be processed elsewhere,+-- while the `flat` transition can be used for more flexible purposes,+-- such as working with samples in memory.+--+-- See <http://msp.org/camcos/2010/5-1/camcos-v5-n1-p04-p.pdf> for the definitive+-- reference of the algorithm.++module Numeric.MCMC.Flat (+    mcmc+  , flat+  , Particle+  , Ensemble+  , Chain++  , module Sampling.Types+  , MWC.create+  , MWC.createSystemRandom+  , MWC.withSystemRandom+  , MWC.asGenIO+  ) where++import Control.Monad (replicateM)+import Control.Monad.Par (NFData)+import Control.Monad.Par.Scheds.Direct hiding (put, get)+import Control.Monad.Par.Combinator (parMap)+import Control.Monad.Primitive (PrimMonad, PrimState, RealWorld)+import Control.Monad.Trans.State.Strict (get, put, execStateT)+import Data.Sampling.Types as Sampling.Types hiding (Chain(..))+import Data.Vector (Vector)+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U+import Pipes (Producer, lift, yield, runEffect, (>->))+import qualified Pipes.Prelude as Pipes+import System.Random.MWC.Probability as MWC++data Chain = Chain {+    chainTarget   :: Target Particle+  , chainPosition :: !Ensemble+  }++instance Show Chain where+  show Chain {..} =+      init+    . filter (`notElem` "[]")+    . unlines+    . V.toList+    . V.map show+    $ chainPosition++type Particle = Vector Double++type Ensemble = Vector Particle++symmetric :: PrimMonad m => Prob m Double+symmetric = fmap transform uniform where+  transform z = 0.5 * (z + 1) ^ (2 :: Int)++stretch :: Particle -> Particle -> Double -> Particle+stretch p0 p1 z = V.zipWith (+) (V.map (* z) p0) (V.map (* (1 - z)) p1)++acceptProb :: Target Particle -> Particle -> Particle -> Double -> Double+acceptProb target particle proposal z =+    lTarget target proposal+  - lTarget target particle+  + log z * (fromIntegral (V.length particle) - 1)++move :: Target Particle -> Particle -> Particle -> Double -> Double -> Particle+move target p0 p1 z zc =+  let proposal = stretch p0 p1 z+      pAccept  = acceptProb target p0 proposal z+  in  if   zc <= min 1 (exp pAccept)+      then proposal+      else p0++execute+  :: PrimMonad m+  => Target Particle+  -> Ensemble+  -> Ensemble+  -> Int+  -> Prob m Ensemble+execute target e0 e1 n = do+  zs  <- replicateM n symmetric+  zcs <- replicateM n uniform+  vjs <- replicateM n (uniformR (1, n))++  let granularity = truncate (fromIntegral n / 2)++      js      = U.fromList vjs+      w0 k    = e0 `V.unsafeIndex` pred k+      w1 k ks = e1 `V.unsafeIndex` pred (ks `U.unsafeIndex` pred k)+++      worker (k, z, zc) = move target (w0 k) (w1 k js) z zc+      result = runPar $+        parMapChunk granularity worker (zip3 [1..n] zs zcs)++  return $ V.fromList result++-- | The 'flat' transition operator for driving a Markov chain over a space+--   of ensembles.+flat+  :: PrimMonad m+  => Transition m Chain+flat = do+  Chain {..} <- get+  let size = V.length chainPosition+      n    = truncate (fromIntegral size / 2)+      e0   = V.slice 0 n chainPosition+      e1   = V.slice n n chainPosition+  result0 <- lift (execute chainTarget e0 e1 n)+  result1 <- lift (execute chainTarget e1 result0 n)+  let ensemble = V.concat [result0, result1]+  put (Chain chainTarget ensemble)++chain :: PrimMonad m => Chain -> Gen (PrimState m) -> Producer Chain m ()+chain = loop where+  loop state prng = do+    next <- lift (MWC.sample (execStateT flat state) prng)+    yield next+    loop next prng++-- | Trace 'n' iterations of a Markov chain and stream them to stdout.+--+--   Note that the Markov chain is defined over the space of ensembles, so+--   you'll need to provide an ensemble of particles for the start location.+--+-- >>> import Numeric.MCMC.Flat+-- >>> import Data.Vector (Vector, toList, fromList)+-- >>> :{+-- >>>   let rosenbrock xs = negate (5  *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)+--             where [x0, x1] = toList xs+-- >>> :}+-- >>> :{+-- >>> let ensemble = fromList [+-- >>>       fromList [negate 1.0, negate 1.0]+-- >>>     , fromList [negate 1.0, 1.0]+-- >>>     , fromList [1.0, negate 1.0]+-- >>>     , fromList [1.0, 1.0]+-- >>>     ]+-- >>> :}+-- >>> withSystemRandom . asGenIO $ mcmc 2 ensemble rosenbrock+-- -1.0,-1.0+-- -1.0,1.0+-- 1.0,-1.0+-- 0.7049046915549257,0.7049046915549257+-- -0.843493377618159,-0.843493377618159+-- -1.1655594505975082,1.1655594505975082+-- 0.5466534497342876,-0.9615123448709006+-- 0.7049046915549257,0.7049046915549257+mcmc :: Int -> Ensemble -> (Particle -> Double) -> Gen RealWorld -> IO ()+mcmc n chainPosition target gen = runEffect $+        chain Chain {..} gen+    >-> Pipes.take n+    >-> Pipes.mapM_ print+  where+    chainTarget = Target target Nothing++-- A parallel map with the specified granularity.+parMapChunk :: NFData b => Int -> (a -> b) -> [a] -> Par [b]+parMapChunk n f xs = concat <$> parMap (map f) (chunk n xs) where+  chunk _ [] = []+  chunk m ys =+    let (as, bs) = splitAt m ys+    in  as : chunk m bs+

@@ -0,0 +1,22 @@+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Main where++import Numeric.MCMC.Flat+import Data.Vector (Vector, toList, fromList)++bnn :: Vector Double -> Double+bnn xs = -0.5 * (x0 ^ 2 * x1 ^ 2 + x0 ^ 2 + x1 ^ 2 - 8 * x0 - 8 * x1) where+  [x0, x1] = toList xs++ensemble :: Ensemble+ensemble = fromList [+    fromList [negate 1.0, negate 1.0]+  , fromList [negate 1.0, 1.0]+  , fromList [1.0, negate 1.0]+  , fromList [1.0, 1.0]+  ]++main :: IO ()+main = withSystemRandom . asGenIO $ mcmc 100 ensemble bnn+

@@ -0,0 +1,22 @@+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Main where++import Numeric.MCMC.Flat+import Data.Vector (Vector, toList, fromList)++rosenbrock :: Vector Double -> Double+rosenbrock xs = negate (5  *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2) where+  [x0, x1] = toList xs++ensemble :: Ensemble+ensemble = fromList [+    fromList [negate 1.0, negate 1.0]+  , fromList [negate 1.0, 1.0]+  , fromList [1.0, negate 1.0]+  , fromList [1.0, 1.0]+  ]++main :: IO ()+main = withSystemRandom . asGenIO $ mcmc 100 ensemble rosenbrock+