diff --git a/flat-mcmc.cabal b/flat-mcmc.cabal --- a/flat-mcmc.cabal +++ b/flat-mcmc.cabal @@ -1,5 +1,5 @@ name: flat-mcmc -version: 1.3.0 +version: 1.4.0 synopsis: Painless general-purpose sampling. homepage: https://github.com/jtobin/flat-mcmc license: MIT @@ -8,7 +8,7 @@ maintainer: jared@jtobin.ca category: Math build-type: Simple -cabal-version: >=1.10 +cabal-version: >= 1.18 description: flat-mcmc is a Haskell library for painless, efficient, general-purpose sampling from continuous distributions. @@ -25,23 +25,23 @@ as a 'flat' transition operator that can be used more generally. . > import Numeric.MCMC.Flat - > import qualified Data.Vector.Unboxed as U (Vector, toList, fromList) - > import qualified Data.Vector as V (fromList) + > import qualified Data.Vector.Unboxed as U (unsafeIndex) > > rosenbrock :: Particle -> Double > rosenbrock xs = negate (5 * (x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2) where - > [x0, x1] = U.toList xs + > x0 = U.unsafeIndex xs 0 + > x1 = U.unsafeIndex xs 1 > - > ensemble :: Ensemble - > ensemble = V.fromList [ - > U.fromList [negate 1.0, negate 1.0] - > , U.fromList [negate 1.0, 1.0] - > , U.fromList [1.0, negate 1.0] - > , U.fromList [1.0, 1.0] + > origin :: Ensemble + > origin = ensemble [ + > particle [negate 1.0, negate 1.0] + > , particle [negate 1.0, 1.0] + > , particle [1.0, negate 1.0] + > , particle [1.0, 1.0] > ] > > main :: IO () - > main = withSystemRandom . asGenIO $ mcmc 12500 ensemble rosenbrock + > main = withSystemRandom . asGenIO $ mcmc 12500 origin rosenbrock Source-repository head Type: git @@ -49,9 +49,16 @@ library default-language: Haskell2010 - ghc-options: -Wall hs-source-dirs: lib - exposed-modules: Numeric.MCMC.Flat + ghc-options: + -Wall + + exposed-modules: + Numeric.MCMC.Flat + + other-modules: + Data.Vector.Extended + build-depends: base > 4 && < 6 , formatting >= 6 && < 7 diff --git a/lib/Data/Vector/Extended.hs b/lib/Data/Vector/Extended.hs new file mode 100644 --- /dev/null +++ b/lib/Data/Vector/Extended.hs @@ -0,0 +1,29 @@ +-- | +-- Module: Data.Vector.Extended +-- Copyright: (c) 2016 Jared Tobin +-- License: MIT +-- +-- Maintainer: Jared Tobin +-- Stability: unstable +-- Portability: ghc + +module Data.Vector.Extended ( + ensemble + , particle + ) where + +import qualified Data.Vector as V (fromList, Vector) +import qualified Data.Vector.Unboxed as U (fromList, Vector) + +-- | A type-specialized alias for Data.Vector.fromList. +-- +-- Use this to create ensembles from lists of particles. +ensemble :: [U.Vector Double] -> V.Vector (U.Vector Double) +ensemble = V.fromList + +-- | A type-specialized alias for Data.Vector.Unboxed.fromList +-- +-- Use this to create particles from lists of doubles. +particle :: [Double] -> U.Vector Double +particle = U.fromList + diff --git a/lib/Numeric/MCMC/Flat.hs b/lib/Numeric/MCMC/Flat.hs --- a/lib/Numeric/MCMC/Flat.hs +++ b/lib/Numeric/MCMC/Flat.hs @@ -36,6 +36,9 @@ , MWC.createSystemRandom , MWC.withSystemRandom , MWC.asGenIO + + , VE.ensemble + , VE.particle ) where import Control.Monad (replicateM) @@ -50,6 +53,7 @@ import qualified Data.Text.IO as T (putStrLn) import Data.Vector (Vector) import qualified Data.Vector as V +import qualified Data.Vector.Extended as VE (ensemble, particle) import qualified Data.Vector.Unboxed as U import Formatting ((%)) import qualified Formatting as F @@ -83,8 +87,20 @@ glue a b = a <> "\n" <> renderParticle b {-# INLINE renderEnsemble #-} +-- | A particle is an n-dimensional point in Euclidean space. +-- +-- You can create a particle by using the 'particle' helper function, or just +-- use Data.Vector.Unboxed.fromList. type Particle = U.Vector Double +-- | An ensemble is a collection of particles. +-- +-- The Markov chain we're interested in will run over the space of ensembles, +-- so you'll want to build an ensemble out of a reasonable number of +-- particles to kick off the chain. +-- +-- You can create an ensemble by using the 'ensemble' helper function, or just +-- use Data.Vector.fromList. type Ensemble = Vector Particle symmetric :: PrimMonad m => Prob m Double @@ -130,7 +146,6 @@ 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) @@ -169,20 +184,20 @@ -- you'll need to provide an ensemble of particles for the start location. -- -- >>> import Numeric.MCMC.Flat --- >>> import Data.Vector (Vector, toList, fromList) +-- >>> import Data.Vector.Unboxed (toList) -- >>> :{ -- >>> 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] +-- >>> let origin = ensemble [ +-- >>> particle [negate 1.0, negate 1.0] +-- >>> , particle [negate 1.0, 1.0] +-- >>> , particle [1.0, negate 1.0] +-- >>> , particle [1.0, 1.0] -- >>> ] -- >>> :} --- >>> withSystemRandom . asGenIO $ mcmc 2 ensemble rosenbrock +-- >>> withSystemRandom . asGenIO $ mcmc 2 origin rosenbrock -- -1.0,-1.0 -- -1.0,1.0 -- 1.0,-1.0 diff --git a/src/Main.hs b/src/Main.hs --- a/src/Main.hs +++ b/src/Main.hs @@ -3,21 +3,20 @@ module Main where import Numeric.MCMC.Flat -import qualified Data.Vector.Unboxed as U (toList, fromList) -import qualified Data.Vector as V (fromList) +import qualified Data.Vector.Unboxed as U (toList) bnn :: Particle -> Double bnn xs = -0.5 * (x0 ^ 2 * x1 ^ 2 + x0 ^ 2 + x1 ^ 2 - 8 * x0 - 8 * x1) where [x0, x1] = U.toList xs -ensemble :: Ensemble -ensemble = V.fromList [ - U.fromList [negate 1.0, negate 1.0] - , U.fromList [negate 1.0, 1.0] - , U.fromList [1.0, negate 1.0] - , U.fromList [1.0, 1.0] +origin :: Ensemble +origin = ensemble [ + particle [negate 1.0, negate 1.0] + , particle [negate 1.0, 1.0] + , particle [1.0, negate 1.0] + , particle [1.0, 1.0] ] main :: IO () -main = withSystemRandom . asGenIO $ mcmc 10000 ensemble bnn +main = withSystemRandom . asGenIO $ mcmc 10000 origin bnn diff --git a/test/BNN.hs b/test/BNN.hs --- a/test/BNN.hs +++ b/test/BNN.hs @@ -3,21 +3,20 @@ module Main where import Numeric.MCMC.Flat -import qualified Data.Vector.Unboxed as U (Vector, toList, fromList) -import qualified Data.Vector as V (fromList) +import qualified Data.Vector.Unboxed as U (toList) bnn :: Particle -> Double bnn xs = -0.5 * (x0 ^ 2 * x1 ^ 2 + x0 ^ 2 + x1 ^ 2 - 8 * x0 - 8 * x1) where [x0, x1] = U.toList xs -ensemble :: Ensemble -ensemble = V.fromList [ - U.fromList [negate 1.0, negate 1.0] - , U.fromList [negate 1.0, 1.0] - , U.fromList [1.0, negate 1.0] - , U.fromList [1.0, 1.0] +origin :: Ensemble +origin = ensemble [ + particle [negate 1.0, negate 1.0] + , particle [negate 1.0, 1.0] + , particle [1.0, negate 1.0] + , particle [1.0, 1.0] ] main :: IO () -main = withSystemRandom . asGenIO $ mcmc 100 ensemble bnn +main = withSystemRandom . asGenIO $ mcmc 100 origin bnn diff --git a/test/Rosenbrock.hs b/test/Rosenbrock.hs --- a/test/Rosenbrock.hs +++ b/test/Rosenbrock.hs @@ -3,21 +3,20 @@ module Main where import Numeric.MCMC.Flat -import qualified Data.Vector.Unboxed as U (Vector, toList, fromList) -import qualified Data.Vector as V (fromList) +import qualified Data.Vector.Unboxed as U (toList) rosenbrock :: Particle -> Double rosenbrock xs = negate (5 *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2) where [x0, x1] = U.toList xs -ensemble :: Ensemble -ensemble = V.fromList [ - U.fromList [negate 1.0, negate 1.0] - , U.fromList [negate 1.0, 1.0] - , U.fromList [1.0, negate 1.0] - , U.fromList [1.0, 1.0] +origin :: Ensemble +origin = ensemble [ + particle [negate 1.0, negate 1.0] + , particle [negate 1.0, 1.0] + , particle [1.0, negate 1.0] + , particle [1.0, 1.0] ] main :: IO () -main = withSystemRandom . asGenIO $ mcmc 100 ensemble rosenbrock +main = withSystemRandom . asGenIO $ mcmc 100 origin rosenbrock