commit 1332db977dc1c0d6205f2f8903d26929a4454b4a
parent 57a93549fd7799d53f59108276bf01b7a28e53de
Author: Jared Tobin <jared@jtobin.ca>
Date:   Sun,  3 Apr 2016 20:34:29 +0700
Merge pull request #1 from jtobin/jt-cleanup
Misc cleanups.
Diffstat:
12 files changed, 352 insertions(+), 217 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -2,7 +2,6 @@
 *swp
 .DS_Store
 dist
-Setup.hs
 *.hi
 *.o
 data
@@ -16,3 +15,8 @@ demos/BNN_Flat
 demos/Rosenbrock_Flat
 demos/Himmelblau_Flat
 demos/SPDE_Flat
+sandbox
+.cabal-sandbox
+cabal.sandbox.config
+.stack-work
+debug
diff --git a/.travis.yml b/.travis.yml
@@ -1 +1,22 @@
-language: haskell
+sudo: false
+language: c
+
+addons:
+  apt:
+    packages:
+    - libgmp-dev
+
+env:
+  STACK_YAML: stack-travis.yaml
+
+before_install:
+- mkdir -p ~/.local/bin
+- export PATH=$HOME/.local/bin:$PATH
+- travis_retry curl -L https://www.stackage.org/stack/linux-x86_64 | tar xz --wildcards --strip-components=1 -C ~/.local/bin '*/stack'
+
+script: stack --no-terminal --install-ghc test --haddock
+
+cache:
+  directories:
+  - $HOME/.stack
+
diff --git a/LICENSE b/LICENSE
@@ -1,29 +1,19 @@
-Copyright (C) 2012 Jared Tobin
+Copyright (c) 2012-2016 Jared Tobin
 
-Redistribution and use in source and binary forms, with or without
-modification, are permitted provided that the following conditions are
-met:
+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:
 
-   1. Redistributions of source code must retain the above copyright
-      notice, this list of conditions, and the following disclaimer.
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
 
-   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.
-
-   3. The name of the author may not be used to endorse or promote
-      products derived from this software without specific prior
-      written permission.
-
-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.
diff --git a/Numeric/MCMC/Flat.hs b/Numeric/MCMC/Flat.hs
@@ -1,169 +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 { _size      :: {-# UNPACK #-} !Int 
-                       , _nEpochs   :: {-# UNPACK #-} !Int  
-                       , _burnIn    :: {-# UNPACK #-} !Int
-                       , _thinEvery :: {-# 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)
-             => ([Double] -> Double)                -- Target to sample
-             -> 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 t e0 e1 n g = do
-    Options _ _ _ _ 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)
-               => ([Double] -> Double)          -- Target to sample
-               -> MarkovChain                   -- State of the Markov chain
-               -> Gen (PrimState m)             -- MWC PRNG
-               -> ViewsOptions m MarkovChain    -- Updated sub-ensemble
-metropolisStep t 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 t e0 e1            n0 g
-    result1 <- executeMoves t e1 (fst result0) n0 g
-
-    return $! 
-      MarkovChain (V.concat $ map fst [result0, result1]) 
-                  (nacc + snd result0 + snd result1)
-{-# INLINE metropolisStep #-}
-
--- | Diffuse through states.
-runChain :: ([Double] -> Double) -- ^ Target to sample
-         -> Options              -- ^ Options of the Markov chain
-         -> MarkovChain          -- ^ Initial state of the Markov chain
-         -> Gen RealWorld        -- ^ MWC PRNG
-         -> IO MarkovChain       -- ^ End state of the Markov chain, wrapped in IO
-runChain target opts initState g 
-    | l == 0 
-        = error "runChain: ensemble must contain at least one particle"
-    | l < (length . V.head) (ensemble initState)
-        = do hPutStrLn stderr $ "runChain: ensemble should be twice as large as "
-                             ++ "the target's dimension.  Continuing anyway."
-             go opts nepochs thinEvery initState g
-    | burnIn < 0 || thinEvery < 0 = error "runChain: nonsensical burn-in or thinning input."
-    | otherwise = go opts nepochs thinEvery initState g
-  where 
-    Options l nepochs burnIn thinEvery _ = opts
-    go o n t !c g0 | n == 0 = hPutStrLn stderr  
-                                 (let nAcc  = accepts c
-                                      total = nepochs * l * length (V.head $ ensemble c)
-                                  in  show nAcc ++ " / " ++ show total ++ " (" ++
-                                      show ((fromIntegral nAcc / fromIntegral total) :: Float) ++ 
-                                      ") proposals accepted") >> return c
-                   | n > (nepochs - burnIn) = do
-                       r <- runReaderT (metropolisStep target c g0) o
-                       go o (n - 1) t r g0 
-                   | n `rem` t /= 0 = do
-                       r <- runReaderT (metropolisStep target c g0) o
-                       go o (n - 1) t r g0
-                   | otherwise = do
-                       r <- runReaderT (metropolisStep target 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 #-}
-
diff --git a/README.md b/README.md
@@ -1,6 +1,43 @@
-# flat-mcmc [](http://travis-ci.org/jtobin/flat-mcmc) 
+# flat-mcmc [](http://travis-ci.org/jtobin/flat-mcmc)
 
-Painless general-purpose sampling.
+*flat-mcmc* is a Haskell library for painless, efficient, general-purpose
+sampling from continuous distributions.
 
-See the *Examples* folder for example usage.
+*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.  Check out
+the paper describing the algorithm
+[here](http://msp.org/camcos/2010/5-1/camcos-v5-n1-p04-p.pdf), and a paper on
+some potential limitations [here](http://arxiv.org/abs/1509.02230).  There is
+also also a robust Python implementation
+[here](http://dan.iel.fm/emcee/current/).
+
+*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.
+
+``` haskell
+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 25000 ensemble rosenbrock
+```
+
+
 
diff --git a/Setup.hs b/Setup.hs
@@ -0,0 +1,7 @@
+module Main (main) where
+
+import Distribution.Simple
+
+main :: IO ()
+main = defaultMain
+
diff --git a/flat-mcmc.cabal b/flat-mcmc.cabal
@@ -1,31 +1,76 @@
--- Initial flat_mcmc.cabal generated by cabal init.  For further 
--- documentation, see http://haskell.org/cabal/users-guide/
-
 name:                flat-mcmc
-version:             0.2.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 25000 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.3, 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
 
diff --git a/lib/Numeric/MCMC/Flat.hs b/lib/Numeric/MCMC/Flat.hs
@@ -0,0 +1,128 @@
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+{-# LANGUAGE RecordWildCards #-}
+
+module Numeric.MCMC.Flat (
+    mcmc
+  , flat
+  , Particle
+  , Ensemble
+
+  , module Sampling.Types
+  , Chain
+  , 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 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 2 worker (zip3 [1..n] zs zcs) -- FIXME granularity option
+
+  return $ V.fromList result
+
+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
+
+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
+
+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
+
diff --git a/stack-travis.yaml b/stack-travis.yaml
@@ -0,0 +1,9 @@
+flags: {}
+packages:
+- '.'
+extra-deps: []
+resolver: lts-5.2
+compiler: ghc-7.10.3
+system-ghc: false
+install-ghc: true
+
diff --git a/stack.yaml b/stack.yaml
@@ -0,0 +1,5 @@
+flags: {}
+packages:
+  - '.'
+extra-deps: []
+resolver: lts-5.2
diff --git a/test/Rosenbrock.hs b/test/Rosenbrock.hs
@@ -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
+
diff --git a/test/Tests.hs b/test/Tests.hs
@@ -0,0 +1,36 @@
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Main where
+
+import Control.Monad
+import Control.Monad.State.Strict
+import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as U
+import Numeric.MCMC.Flat
+
+lRosenbrock :: Density
+lRosenbrock xs =
+  let [x0, x1] = U.toList xs
+  in  (-1) * (5 * (x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
+
+defaultEnsemble :: Ensemble
+defaultEnsemble = V.fromList $ map U.fromList 
+  [[0.1, 0.5], [0.8, 0.1], [1.0, 0.2], [0.9, 0.8], [-0.2, 0.3], [-0.1, 0.9]]
+
+opts :: Options
+opts = Options 10
+
+origin :: Chain
+origin = Chain {
+    logObjective = lRosenbrock
+  , ensemble     = defaultEnsemble
+  , iterations   = 0
+  , accepts      = 0
+  } 
+
+-- cabal test --show-details=streaming
+main :: IO ()
+main = withSystemRandom . asGenIO $ \g -> do
+  trace <- sample (replicateM 5000 (flatGranular opts) `evalStateT` origin) g
+  print trace
+