diff --git a/ChangeLog.md b/ChangeLog.md
new file mode 100644
--- /dev/null
+++ b/ChangeLog.md
@@ -0,0 +1,3 @@
+# Changelog for markov
+
+## Unreleased changes
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Author name here (c) 2019
+
+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 Author name here 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,66 @@
+# Markov Tutorial
+
+Let Xn denote the nth state of a Markov chain with state space ℕ.
+For x ≠ 0 define transition probabilities
+
+p(x,0) = q,
+
+p(x,x) = r, and
+
+p(x,x+1) = s.
+
+When x = 0, let
+p(x,0) = q+r,
+p(x,x+1) = s.
+Let p(x,y) = 0 in all other cases.
+Suppose we wanted to find
+P\[Xn = j ∩ d = k],
+where d denotes the number of transitions from a positive integer to zero.
+There are three values we need to track —
+extinctions, probability, and state.
+Extinctions add a value to a counter each time they happen
+and the counter takes integral values,
+so they can be represented by `Sum Int`.
+Probabilities are multiplied each step,
+and added when duplicate steps are combined.
+We want decimal probabilities, so
+we can represent this with `Product Rational`.
+We will make a new type for the state.
+
+```haskell
+data Extinction = Extinction Int
+ deriving Generic
+ deriving newtype (Eq, Num, Show)
+ deriving anyclass Grouping
+```
+
+All that remains is to make an instance of `Markov`.
+
+```haskell
+instance Markov (Sum Int, Product Rational) Extinction where
+ transition x = case state x of
+ 0 -> [ 0 >*< (q+r) >*< id
+ , 0 >*< s >*< (+1) ]
+ _ -> [ 1 >*< q >*< const 0
+ , 0 >*< r >*< id
+ , 0 >*< s >*< (+1) ]
+ where q = 0.1; r = 0.3; s = 0.6
+```
+
+We can now easily see a list of states, deaths, and the probabilities.
+
+
+
+> chain [pure 0 :: Sum Int :* Product Rational :* Extinction] !! 3
+
+((0,8 % 125),0)
+((0,111 % 500),1)
+((1,51 % 500),0)
+((0,9 % 25),2)
+((1,9 % 250),1)
+((0,27 % 125),3
+
+
+This means that starting from a state of zero,
+after three time steps there is a 51/500 chance
+that the state is zero and there has been one death.
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/markov-realization.cabal b/markov-realization.cabal
new file mode 100644
--- /dev/null
+++ b/markov-realization.cabal
@@ -0,0 +1,38 @@
+cabal-version: 1.12
+name: markov-realization
+version: 0.1.0
+license: BSD3
+license-file: LICENSE
+copyright: 2019 Alex Loomis
+maintainer: atloomis@math.arizona.edu
+author: Alex Loomis
+homepage: https://github.com/alexloomis/markov
+bug-reports: https://github.com/alexloomis/markov/issues
+synopsis: Realizations of Markov chains.
+description:
+ Please see the README on GitHub at
+category: Statistics
+build-type: Simple
+extra-source-files:
+ README.md
+ ChangeLog.md
+
+source-repository head
+ type: git
+ location: git://github.com/alexloomis/markov.git
+
+library
+ exposed-modules:
+ Markov
+ Markov.Examples
+ Markov.Instances
+ hs-source-dirs: src
+ other-modules:
+ Paths_markov_realization
+ default-language: Haskell2010
+ build-depends:
+ base >=4.7 && <5,
+ contravariant >=1.5.1 && <1.6,
+ discrimination ==0.4.*,
+ generic-deriving >=1.12.4 && <1.13,
+ MonadRandom >=0.5.1.1 && <0.6
diff --git a/src/Markov.hs b/src/Markov.hs
new file mode 100644
--- /dev/null
+++ b/src/Markov.hs
@@ -0,0 +1,219 @@
+{-# LANGUAGE MultiParamTypeClasses, GeneralizedNewtypeDeriving, DeriveGeneric,
+DeriveAnyClass, DerivingStrategies, TypeOperators #-}
+{-|
+Module : Markov
+Description : Realization of Markov processes with known parameters.
+Maintainer : atloomis@math.arizona.edu
+Stability : experimental
+
+Three type classes for deterministically analyzing
+Markov chains with known parameters.
+'Markov0' is intended to list possible outcomes,
+'Markov' should allow for more sophisticated analysis,
+and 'MultiMarkov' is intended to make implementing
+hidden Markov models easier.
+See "Examples" for examples.
+See README for a detailed description.
+-}
+module Markov (
+ -- *Markov0
+ Markov0 (..)
+ -- *Markov
+ , Markov (..)
+ -- *MultiMarkov
+ , randomProduct
+ , randomPath
+ , MultiMarkov (..)
+ -- *Combine
+ , Combine (..)
+ , Merge (..)
+ , Sum (..)
+ , Product (..)
+ -- *Misc
+ , (:*)
+ , (>*<)
+ , fromLists
+ -- *Testing
+ ) where
+
+import Markov.Instances ()
+import Control.Applicative ((<**>))
+import Generics.Deriving (Generic)
+import Data.Discrimination (Grouping, grouping)
+import qualified Data.Discrimination as DD
+import qualified Data.List as DL
+import qualified Data.List.NonEmpty as NE
+import qualified Data.Functor.Contravariant as FC
+import qualified Control.Monad.Random as MR
+
+---------------------------------------------------------------
+-- Markov0
+---------------------------------------------------------------
+
+-- |A basic implementation of Markov chains.
+class (Eq m) => Markov0 m where
+ -- |The transition functions from a state.
+ transition0 :: m -> [m -> m]
+ step0 :: m -> [m]
+ -- |Iterated steps.
+ chain0 :: [m] -> [[m]]
+ step0 x = fmap ($ x) (transition0 x)
+ chain0 = DL.iterate' $ DL.nub . concatMap step0
+
+---------------------------------------------------------------------------------------
+-- Markov
+---------------------------------------------------------------------------------------
+
+-- |An implementation of Markov chains.
+-- To speed up @chain@, try instead:
+--
+-- > chain = DL.iterate' $ map summarize' . NE.group . DL.sort . concatMap step
+-- > where summarize' xs@((_,b)NE.:|_) = (summarize . fmap fst $ xs, b)
+class (Combine t, Grouping t, Grouping m, Monoid t) => Markov t m where
+ transition :: m -> [(t, m -> m)]
+ step :: (t,m) -> [(t,m)]
+ chain :: [(t,m)] -> [[(t,m)]]
+ step x = fmap (x <**>) (transition $ snd x)
+ -- |Iterated steps, with equal states combined using 'summarize' operation.
+ chain = DL.iterate' $ map (summarize' . NE.fromList)
+ . DD.group . concatMap step
+ where summarize' xs@((_,b)NE.:|_) = (summarize . fmap fst $ xs, b)
+ -- WARNING: DD.group does not currently respect equivalence classes.
+
+---------------------------------------------------------------------------------------
+-- Multi-Transition Markov
+---------------------------------------------------------------------------------------
+
+-- |An implementation of Markov chains that allows multi-transition steps.
+class (Combine m, Grouping m, Semigroup m) => MultiMarkov m where
+ multiTransition :: m -> [m -> [m]]
+ multiStep :: m -> [m]
+ multiChain :: [m] -> [[m]]
+ multiStep x = foldr phi [x] (multiTransition x)
+ where phi f = concatMap (delta f)
+ delta f y = map (y <>) (f y)
+ multiChain = DL.iterate' $ map (summarize . NE.fromList)
+ . DD.group . concatMap multiStep
+
+---------------------------------------------------------------------------------------
+-- Combine
+---------------------------------------------------------------------------------------
+
+-- |Within equivalence classes, @combine@ should be associative,
+-- commutative, and should be idempotent up to equivalence.
+-- I.e. if @x == y == z@,
+--
+-- prop> (x `combine` y) `combine` z = x `combine` (y `combine` z)
+-- prop> x `combine` y = y `combine` x
+-- prop> x `combine` x == x
+class Combine a where
+ combine :: a -> a -> a
+ summarize :: NE.NonEmpty a -> a
+ summarize (a NE.:| b) = foldr combine a b
+
+instance (Combine a, Combine b) => Combine (a,b) where
+ combine (w,x) (y,z) = (combine w y, combine x z)
+
+instance (Combine a, Combine b, Combine c) => Combine (a,b,c) where
+ combine (a,w,x) (b,y,z) = (combine a b, combine w y, combine x z)
+
+---------------------------------------------------------------------------------------
+-- Easier way to write nested 2-tuples
+---------------------------------------------------------------------------------------
+
+-- |Easier way to write nested 2-tuples.
+type a :* b = (a,b)
+-- |Easier way to write nested 2-tuples,
+-- since @a >*\< b >*\< c >*< d@
+-- is much easier to read than
+-- @(((a,b),c),d)@.
+-- Left associative, binds weaker than @+@
+-- but stronger than @==@.
+(>*<) :: a -> b -> a :* b
+a >*< b = (a,b)
+infixl 5 >*<
+
+---------------------------------------------------------------------------------------
+-- Merge
+---------------------------------------------------------------------------------------
+
+-- Does not group to combine unless equal.
+-- |Values from a 'Monoid' which have the respective
+-- binary operation applied each step.
+-- E.g., strings with concatenation.
+newtype Merge a = Merge a
+ deriving (Eq, Generic)
+ deriving newtype (Semigroup, Monoid, Enum, Num, Fractional, Show)
+ deriving anyclass Grouping
+
+instance Combine (Merge a) where
+ combine = const
+
+---------------------------------------------------------------------------------------
+-- Sum
+---------------------------------------------------------------------------------------
+
+-- |Values which are added each step.
+-- E.g., number of times a red ball is picked from an urn.
+newtype Sum a = Sum a
+ deriving Generic
+ deriving newtype (Eq, Enum, Num, Fractional, Show)
+ deriving anyclass Grouping
+
+instance Combine (Sum a) where
+ combine = const
+
+instance Num a => Semigroup (Sum a) where
+ x <> y = x + y
+
+instance Num a => Monoid (Sum a) where
+ mempty = 0
+
+---------------------------------------------------------------------------------------
+-- Product
+---------------------------------------------------------------------------------------
+
+-- Does not effect equality of tuple,
+-- @combine x y = x + y@.
+-- |Values which are multiplied each step,
+-- and combined additively for equal states.
+-- E.g., probabilities.
+newtype Product a = Product a
+ deriving Generic
+ deriving newtype (Num, Fractional, Enum, Show)
+
+instance Grouping (Product a) where
+ grouping = FC.contramap (const ()) grouping
+
+-- This causes Data.List.group to act more like Data.Discrimination.group
+instance Eq (Product a) where
+ _ == _ = True
+
+instance Num a => Combine (Product a) where
+ combine = (+)
+
+instance Num a => Semigroup (Product a) where
+ x <> y = x * y
+
+instance Num a => Monoid (Product a) where
+ mempty = 1
+
+---------------------------------------------------------------------------------------
+-- Misc
+---------------------------------------------------------------------------------------
+
+-- |Randomly choose from a list by probability.
+randomProduct :: (Real a, MR.MonadRandom m) => [(a, b)] -> m (a, b)
+randomProduct xs = MR.fromList . map (\x -> (x, toRational $ fst x)) $ xs
+
+-- |Returns a single realization of a Markov chain.
+randomPath :: (Markov a b, Real a, MR.RandomGen g) => (a,b) -> g -> [(a,b)]
+randomPath x g = map (flip MR.evalRand g) . iterate (>>= (randomProduct . step)) $ pure x
+
+-- |Create a transition function from a transition matrix.
+-- If [[a]] is an n x n matrix, length [b] should be n.
+fromLists :: Eq b => [[a]] -> [b] -> b -> [(a, b -> b)]
+fromLists matrix states b = case DL.elemIndex b states of
+ Nothing -> []
+ Just n -> zip (matrix!!n) toState
+ where toState = map const states
diff --git a/src/Markov/Examples.hs b/src/Markov/Examples.hs
new file mode 100644
--- /dev/null
+++ b/src/Markov/Examples.hs
@@ -0,0 +1,422 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, DeriveGeneric, DeriveAnyClass,
+DerivingStrategies, GeneralizedNewtypeDeriving, FlexibleInstances, TypeOperators #-}
+{-|
+Module : Examples
+Description : Examples of Markov chains implemented using "Markov".
+Maintainer : atloomis@math.arizona.edu
+Stability : experimental
+
+Several examples of Markov chains.
+It is probably more helpful to read the source code than the Haddock documentation.
+-}
+module Markov.Examples ( FromMatrix (..)
+ , Simple (..)
+ , Urn (..)
+ , Extinction (..)
+ , Tidal (..)
+ , Room (..)
+ , FillBin
+ , initial
+ , expectedLoss
+ ) where
+
+import Markov
+import Generics.Deriving (Generic)
+import Data.Discrimination (Grouping)
+
+---------------------------------------------------------------
+-- From a matrix
+---------------------------------------------------------------
+
+-- |An example defined from a matrix.
+--
+-- >>> chain [pure 't' :: Product Double :* Char] !! 100
+-- [ (0.5060975609756099,'a')
+-- , (0.201219512195122,'t')
+-- , (0.29268292682926833,'l') ]
+newtype FromMatrix = FromMatrix Char
+ deriving Generic
+ deriving newtype (Eq, Show)
+ deriving anyclass Grouping
+
+instance Markov (Product Double) FromMatrix where
+ transition = let mat = [ [0.4, 0.3, 0.3]
+ , [0.2, 0.1, 0.7]
+ , [0.9, 0.1, 0.0] ]
+ chars = map FromMatrix ['a','t','l']
+ in fromLists mat chars
+
+---------------------------------------------------------------
+-- Simple random walk
+---------------------------------------------------------------
+
+-- |A simple random walk.
+-- Possible outcomes of the first three steps:
+--
+-- >>> take 3 $ chain0 [Simple 0]
+-- [ [0]
+-- , [-1,1]
+-- , [-2,0,2]]
+--
+-- Probability of each outcome:
+--
+-- >>> take 3 $ chain [pure 0 :: Product Double :* Simple]
+-- [ [(1.0,0)]
+-- , [(0.5,-1),(0.5,1)]
+-- , [(0.25,-2),(0.5,0),(0.25,2)] ]
+--
+-- Number of ways to achieve each outcome:
+--
+-- >>> take 3 $ chain [pure 0 :: Product Int :* Simple]
+-- [ [(1,0)]
+-- , [(1,-1),(1,1)]
+-- , [(1,-2),(2,0),(1,2)] ]
+--
+-- Number of times @pred@ was applied,
+-- allowing steps in place (@id@)
+-- for more interesting output:
+--
+-- >>> chain [pure 0 :: Sum Int :* Simple] !! 2
+-- [ (2,-2)
+-- , (1,-1)
+-- , (1,0)
+-- , (0,0)
+-- , (0,1)
+-- , (0,2) ]
+
+newtype Simple = Simple Int
+ deriving Generic
+ deriving newtype (Num, Enum, Eq, Ord, Show)
+ deriving anyclass Grouping
+
+instance Markov0 Simple where
+ transition0 _ = [pred, succ]
+
+instance Markov (Product Double) Simple where
+ transition _ = [ 0.5 >*< pred
+ , 0.5 >*< succ ]
+
+instance Markov (Product Int) Simple where
+ transition _ = [ 1 >*< pred
+ , 1 >*< succ ]
+
+instance Markov (Sum Int) Simple where
+ transition _ = [ 1 >*< pred
+ , 0 >*< id
+ , 0 >*< succ ]
+ -- = [ 1 >*< pred
+ -- , pure id
+ -- , pure succ ]
+
+---------------------------------------------------------------
+-- Urn model
+---------------------------------------------------------------
+
+-- |An urn contains balls of two colors.
+-- At each step, a ball is chosen uniformly at random from the urn
+-- and a ball of the same color is added.
+newtype Urn = Urn (Int,Int)
+ deriving Generic
+ deriving newtype (Eq, Ord, Show)
+ deriving anyclass Grouping
+
+instance Markov (Product Double) Urn where
+ transition x = [ probLeft x >*< addLeft
+ , 1 - probLeft x >*< addRight ]
+
+addLeft :: Urn -> Urn
+addLeft (Urn (a,b)) = Urn (a+1,b)
+
+addRight :: Urn -> Urn
+addRight (Urn (a,b)) = Urn (a,b+1)
+
+probLeft :: Fractional a => Urn -> a
+probLeft (Urn (a,b)) = (fromIntegral a)/(fromIntegral $ a + b)
+
+---------------------------------------------------------------
+-- Tutorial
+---------------------------------------------------------------
+
+-- |This is the chain from the README.
+newtype Extinction = Extinction Int
+ deriving Generic
+ deriving newtype (Eq, Num, Show)
+ deriving anyclass Grouping
+
+instance Markov (Sum Int, Product Rational) Extinction where
+ transition x = case x of
+ 0 -> [ 0 >*< (q+r) >*< id
+ , 0 >*< s >*< (+1) ]
+ _ -> [ 1 >*< q >*< const 0
+ , 0 >*< r >*< id
+ , 0 >*< s >*< (+1) ]
+ where q = 0.1; r = 0.3; s = 0.6
+
+-- This is equivalent to the definition above.
+instance Combine Extinction where
+ combine = const
+
+instance Semigroup Extinction where
+ (<>) = flip const
+
+instance MultiMarkov (Sum Int :* Product Rational :* Extinction) where
+ multiTransition _ = [trans]
+ where trans ((_,_),z) = case z of
+ 0 -> [ 0 >*< (q+r) >*< 0
+ , 0 >*< s >*< 1 ]
+ x -> [ 1 >*< q >*< 0
+ , 0 >*< r >*< x
+ , 0 >*< s >*< x+1 ]
+ where q = 0.1; r = 0.3; s = 0.6
+
+---------------------------------------------------------------
+-- More complex random walk
+---------------------------------------------------------------
+
+-- |A time inhomogenous random walk that vaguely models tides
+-- by periodically switching directions
+-- and falling back from a shore at the origin.
+data Tidal = Tidal { time :: Double
+ , position :: Int }
+ deriving (Eq, Ord, Show, Generic)
+ deriving anyclass Grouping
+
+instance Markov (Product Double) Tidal where
+ transition tw = [ probRight tw >*< stepPos (+1)
+ , 1 - (probRight tw) >*< stepPos (flip (-) 1) ]
+
+stepPos :: (Int -> Int) -> Tidal -> Tidal
+stepPos f tw = Tidal (time tw + 1) (f $ position tw)
+
+probRight :: Tidal -> Product Double
+probRight tw = Product $ timeBias * positionBias
+ where timeBias = (1 + sin (2 * pi * (time tw) / stepsPerCycle))/2
+ positionBias
+ | position tw >= 0 = 1 / steepness
+ | otherwise = 1
+ stepsPerCycle = 10
+ steepness = 1.3 -- Double from 1 (flat) to +infty
+
+---------------------------------------------------------------
+-- Hidden Markov Model
+---------------------------------------------------------------
+
+-- |A hidden Markov model.
+--
+-- >>> filter (\((_,Merge xs),_) -> xs == "aaa") $ multiChain [1 >*< Merge "" >*< 1 :: Product Rational :* Merge String :* Room] !! 3
+-- [ ((3243 % 200000,"aaa"),Room 1)
+-- , ((9729 % 500000,"aaa"),Room 2)
+-- , ((4501 % 250000,"aaa"),Room 3) ]
+--
+-- Given that all three tokens recieved were @"a"@,
+-- there is a probability of approximately @0.34@
+-- that the current room is @Room 3@.
+newtype Room = Room Int
+ deriving (Generic, Show)
+ deriving newtype (Eq, Num)
+ deriving anyclass Grouping
+
+instance Semigroup Room where
+ (<>) = flip const
+
+instance Combine Room where
+ combine = const
+
+-- Note that changeState is applied before giveToken.
+-- In spirit, we have stepj = giveToken . changeState
+instance MultiMarkov (Product Rational :* Merge String :* Room) where
+ multiTransition _ = [giveToken, changeState]
+ where changeState ((_,_),z) = case z of
+ 1 -> [ 0.3 >*< mempty >*< 1
+ , 0.6 >*< mempty >*< 2
+ , 0.1 >*< mempty >*< 3 ]
+ 2 -> [ 1.0 >*< mempty >*< 3 ]
+ 3 -> [ 0.3 >*< mempty >*< 1
+ , 0.6 >*< mempty >*< 2
+ , 0.1 >*< mempty >*< 3 ]
+ _ -> error "State out of bounds in transitionk"
+ giveToken ((_,_),z) = case z of
+ 1 -> [ 0.5 >*< Merge "a" >*< 1
+ , 0.5 >*< Merge "b" >*< 1 ]
+ 2 -> [ 0.3 >*< Merge "a" >*< 2
+ , 0.7 >*< Merge "b" >*< 2 ]
+ 3 -> [ 0.4 >*< Merge "a" >*< 3
+ , 0.4 >*< Merge "b" >*< 3
+ , 0.2 >*< Merge "c" >*< 3 ]
+ _ -> error "State out of bounds in transitionk"
+
+---------------------------------------------------------------
+-- Yet more complex example
+---------------------------------------------------------------
+
+-- |Represents bins with free slots and items.
+type Bin = (Open,Full)
+type Index = Int
+-- |Represents space between bins where they can expand.
+type Gap = Int
+type Full = Int
+type Open = Int
+type Trans = FillBin -> FillBin
+
+-- |A collection of bins with gaps between them.
+-- At each step an empty space is chosen
+-- form a bin or from a gap.
+-- If it is in a bin, the space is filled.
+-- If it is in a gap, it is assigned to an adjacent bin,
+-- which expands to contain it and any intervening spaces,
+-- and then the space filled.
+data FillBin = End Gap | Ext Gap Bin FillBin deriving (Eq, Ord, Generic, Grouping)
+
+instance Show FillBin where
+ show (Ext g b s) = show g ++ " " ++ show b ++ " " ++ show s
+ show (End g) = show g
+
+instance Markov (Product Double) FillBin where
+ transition x = case probId x of
+ 0 -> filter (\(Product y,_) -> y /= 0) -- Careful, Product _ == Product _ = True
+ $ [probAdd i x >*< addItem i | i <- indices]
+ ++ [probGrowL i x >*< addItem i . growLeft j i
+ | i <- indices, j <- [1..gapN (i-1) x]]
+ ++ [probGrowR i x >*< addItem i . growRight j i
+ | i <- indices, j <- [1..gapN i x]]
+ 1 -> [pure id]
+ _ -> error "Pattern not matched in transition"
+ where indices = [1..size x]
+
+-- |>>> fBFromLists [1,3,5,10] [(3,5),(9,9),(8,3)]
+-- 1 (3,5) 3 (9,9) 5 (8,3) 10
+fBFromLists :: [Gap] -> [Bin] -> FillBin
+fBFromLists gaps bins = case (gaps,bins) of
+ (g:_ , [] ) -> End g
+ ([g] , _ ) -> End g
+ (g:gs , b:bs) -> Ext g b $ fBFromLists gs bs
+ ([] , _ ) -> End 0
+
+-- |Create state where all bins start as (0,0).
+--
+-- >>> initial [5,7,0]
+-- 5 (0,0) 7 (0,0) 0
+initial :: [Int] -> FillBin
+initial gs = fBFromLists gs $ repeat (0,0)
+
+-- |The number of bins.
+size :: FillBin -> Int
+size x = case x of
+ End _ -> 0
+ Ext _ _ s -> 1 + size s
+
+-- |The bins of a state.
+getBins :: FillBin -> [Bin]
+getBins x = case x of
+ End _ -> []
+ Ext _ b s -> b:getBins s
+
+-- |The open values of a state.
+getOpen :: FillBin -> [Open]
+getOpen x = map fst $ getBins x
+
+-- |The open value of the Nth bin.
+openN :: Index -> FillBin -> Open
+openN i x = (getOpen x)!!(i-1)
+
+-- |The full values of a state.
+getFull :: FillBin -> [Full]
+getFull x = map snd $ getBins x
+
+-- |The full value of the Nth bin.
+fullN :: Index -> FillBin -> Full
+fullN i x = (getFull x)!!(i-1)
+
+-- |The gap values of a state.
+getGap :: FillBin -> [Gap]
+getGap x = case x of
+ End g -> [g]
+ Ext g _ s -> g:getGap s
+
+-- |Warning! Indexed from zero!
+gapN :: Index -> FillBin -> Gap
+gapN i x = (getGap x)!!i
+
+-- |The command @iApply i f s@ is analagous to
+-- @take i s ++ f (drop i s)@.
+iApply :: Trans -> Index -> Trans
+iApply f idx x = case (idx,x) of
+ (1, y) -> f y
+ (i, Ext g b s) -> Ext g b $ iApply f (i-1) s
+ _ -> error "Pattern not matched in iApply"
+
+-- |Add an item to the ith bin.
+addItem :: Index -> Trans
+addItem = iApply h
+ where h (Ext g (o,f) s) = Ext g (o-1,f+1) s
+ h _ = error "pattern not matched in h in addItem"
+
+-- |Expand the ith bin to the left by j.
+-- The Markov chain will use @addItem i . growLeft j i@.
+growLeft :: Int -> Index -> Trans
+growLeft j = iApply h
+ where h (Ext g (o,f) s) = Ext (g-j) (o+j,f) s
+ h _ = error "pattern not matched in h in growLeft"
+
+growRight :: Int -> Index -> Trans
+growRight j = iApply h
+ where h (Ext g (o,f) s) = Ext g (o+j,f) (shrink s)
+ h _ = error "pattern not matched in h in growRight"
+ shrink s = case s of
+ End g -> End (g-j)
+ Ext g b t -> Ext (g-j) b t
+
+-- |The sum of all open slots in bins and gaps.
+slots :: FillBin -> Int
+slots x = sum $ getGap x ++ getOpen x
+
+-- |The probability that a state returns to itself.
+probId :: Num a => FillBin -> a
+probId x = case slots x == 0 of
+ True -> 1
+ False -> 0
+
+divInt :: (Integral a, Integral b, Fractional c) => a -> b -> c
+divInt x y = (fromIntegral x)/(fromIntegral y)
+
+-- |The probability that the ith bin gains an item.
+probAdd :: Fractional a => Index -> FillBin -> a
+probAdd i x = openN i x `divInt` slots x
+
+-- |The probability that the ith bin expands to the left.
+probGrowL :: Fractional a => Index -> FillBin -> a
+probGrowL i x = case test of
+ True -> 1 `divInt` slots x
+ False -> 0
+ where test = i == 1 || fullN i x < fullN (i-1) x
+
+-- |The probability that the ith bin expands to the right.
+probGrowR :: Fractional a => Index -> FillBin -> a
+probGrowR i x = case test of
+ True -> 1 `divInt` slots x
+ False -> 0
+ where test = i == size x || fullN i x <= fullN (i+1) x
+
+---------------------------------------------------------------
+-- Several functions to help study the previous process
+---------------------------------------------------------------
+
+-- |The \(l^2\) distance between a finished state
+-- and a state with perfectly balanced bins.
+individualLoss :: Fractional a => FillBin -> a
+individualLoss x = sum . map f . getFull $ x
+ where f y = (fromIntegral y - ideal)^2
+ ideal = sum (getFull x) `divInt` size x
+
+probLoss :: Fractional a => (Product a, FillBin) -> a
+probLoss (Product x, y) = x * individualLoss y
+
+-- |Expected loss of a set of pstates of @['FillBin']@.
+-- Loss is the \(l^2\) distance between a finished state
+-- and a state with perfectly balanced bins.
+--
+-- >>> expectedLoss [pure $ initial [1,0,3] :: Product Double :* FillBin]
+-- 2.0
+expectedLoss :: (Fractional a, Markov (Product a) FillBin) => [Product a :* FillBin] -> a
+expectedLoss xs = sum . map probLoss $ (chain xs) !! idx
+ where idx = slots . snd . head $ xs
diff --git a/src/Markov/Instances.hs b/src/Markov/Instances.hs
new file mode 100644
--- /dev/null
+++ b/src/Markov/Instances.hs
@@ -0,0 +1,9 @@
+module Markov.Instances where
+
+import Data.Discrimination (Grouping, grouping)
+import Data.Functor.Contravariant (contramap)
+import Data.Functor.Contravariant.Divisible (conquer)
+import GHC.Float (castFloatToWord32, castDoubleToWord64)
+
+instance Grouping Float where grouping = contramap castFloatToWord32 grouping
+instance Grouping Double where grouping = contramap castDoubleToWord64 grouping