Haskell implementations of Hidden Markov Models & related algorithms from AIMA book
modern nix setup, proper haddock documentations, bug fixes
+225
-115
+4
-1
.gitignore
+4
-1
.gitignore
+7
-11
default.nix
+7
-11
default.nix
···
1
-
{ nixpkgs ? import <nixpkgs> {}, compiler ? "default", withHoogle ? false }:
1
+
{ nixpkgs ? import <nixpkgs> {}, compiler ? "default"}:
2
2
let
3
3
inherit (nixpkgs) pkgs;
4
4
haskellPackages = if compiler == "default"
5
5
then pkgs.haskellPackages
6
6
else pkgs.haskell.packages.${compiler};
7
7
8
-
hPkgs = if withHoogle
9
-
then
10
-
haskellPackages // rec {
11
-
ghc = haskellPackages.ghc // { withPackages = haskellPackages.ghc.withHoogle; };
12
-
ghcWithPackages = ghc.withPackages;
13
-
}
14
-
else
15
-
haskellPackages;
16
-
17
8
in
18
-
hPkgs.callPackage ./hmm.nix {}
9
+
{
10
+
hmm = haskellPackages.developPackage {
11
+
returnShellEnv = false;
12
+
root = ./.;
13
+
};
14
+
}
-12
hmm.nix
-12
hmm.nix
···
1
-
{ mkDerivation, base, combinatorial, hmatrix, linear, mtl, papa
2
-
, stdenv, transformers, vector
3
-
}:
4
-
mkDerivation {
5
-
pname = "hmm";
6
-
version = "0.1.0.0";
7
-
src = ./.;
8
-
libraryHaskellDepends = [
9
-
base combinatorial hmatrix linear mtl papa transformers vector
10
-
];
11
-
license = stdenv.lib.licenses.asl20;
12
-
}
+9
-3
shell.nix
+9
-3
shell.nix
···
1
1
{ nixpkgs ? import <nixpkgs> {}, compiler ? "default", withHoogle ? false}:
2
2
let
3
3
inherit (nixpkgs) pkgs;
4
-
drv = import ./default.nix { inherit nixpkgs compiler withHoogle; };
5
-
drvWithTools = pkgs.haskell.lib.addBuildDepends drv [ pkgs.cabal-install ];
4
+
haskellPackages = if compiler == "default"
5
+
then pkgs.haskellPackages
6
+
else pkgs.haskell.packages.${compiler};
7
+
8
+
project = import ./. { inherit nixpkgs compiler; };
6
9
in
7
-
if pkgs.lib.inNixShell then drvWithTools.env else drvWithTools
10
+
haskellPackages.shellFor {
11
+
withHoogle = true;
12
+
packages = p: [ project.hmm ];
13
+
}
+12
-5
src/cHMM.hs
+12
-5
src/cHMM.hs
···
1
1
{-# LANGUAGE DataKinds #-}
2
2
{-# LANGUAGE TypeOperators #-}
3
3
{-# LANGUAGE RankNTypes #-}
4
-
{-# LANGUAGE TypeApplications #-}
4
+
---------------------------------
5
+
-- |
6
+
-- Module : CHMM
7
+
-- Copyright : (C) Kaushik Chakraborty, 2019
8
+
-- License : Apache v2 (see the file LICENSE)
9
+
-- Maintainer : Kaushik Chakraborty <git@kaushikc.org>
10
+
-- Stability : experimental
11
+
--
12
+
---------------------------------
5
13
6
14
module CHMM where
7
15
8
-
import Prelude(undefined)
9
16
import Papa
10
-
11
17
import GHC.TypeNats
12
18
import Data.Proxy
13
19
import Linear.V
···
17
23
-----------------------
18
24
19
25
trajectory :: forall n t a . (KnownNat n, KnownNat t) => Proxy t -> V n a -> V (n ^ t) (V n a)
20
-
trajectory p = V . Vec.fromList . fmap (V . Vec.fromList) . variateRep t . toList
26
+
trajectory p = V . Vec.fromList . fmap (V . Vec.fromList) . variateRep t' . toList
21
27
where
22
-
t = fromIntegral $ natVal p
28
+
t' = fromIntegral $ natVal p
29
+
+193
-83
src/hmm.hs
+193
-83
src/hmm.hs
···
3
3
{-# LANGUAGE ScopedTypeVariables #-}
4
4
{-# LANGUAGE TemplateHaskell #-}
5
5
{-# LANGUAGE TypeOperators #-}
6
+
---------------------------------
7
+
-- |
8
+
-- Module : HMM
9
+
-- Copyright : (C) Kaushik Chakraborty, 2019
10
+
-- License : Apache v2 (see the file LICENSE)
11
+
-- Maintainer : Kaushik Chakraborty <git@kaushikc.org>
12
+
-- Stability : experimental
13
+
--
14
+
---------------------------------
15
+
6
16
module HMM where
7
17
8
-
-- import Prelude (undefined)
9
-
-- import System.IO (print)
18
+
import Prelude (undefined)
19
+
import System.IO (print)
10
20
import Text.Printf (printf)
11
21
12
22
import Papa
···
25
35
import qualified Numeric.LinearAlgebra as HMat
26
36
-----------------------------------------------------------------------------------
27
37
28
-
-- | Type synonyms
38
+
-- * Type synonyms
29
39
type R = Double
30
40
type TransitionModel (s :: Nat) a = V s (V s a)
31
41
type SensorDiagonal (s :: Nat) a = V s (V s a)
···
33
43
type Message (s :: Nat) a = V s (V 1 a)
34
44
type Distribution (s :: Nat) a = V s a
35
45
36
-
-- | Utility Functions
46
+
-- * HMM
47
+
48
+
-- | Hidden Markov Models
49
+
data HMM (s :: Nat) (t :: Nat) a b = HMM {
50
+
-- | prior distribution which is used as initial forward message
51
+
_prior :: Message s a,
52
+
-- | transition model with @s@ states as @sxs@ matrix
53
+
_tModel :: TransitionModel s a,
54
+
-- | evidence value vector, each index maps to the corresponding evidence values
55
+
_sDist :: V t b,
56
+
-- | sensor model with @t@ evidence values having @sxs@ diagonal matrix capturing their ditributions for each state @s@
57
+
_sModel :: SensorModel t s a
58
+
} deriving (Show)
59
+
60
+
-- |
61
+
-- === __HMM lenses__
62
+
--
63
+
makeLenses ''HMM
64
+
65
+
-- ** Smart constructors
66
+
mkHMM_ :: (KnownNat s, KnownNat t) => Maybe (Message s R) -- ^ Prior distribution on the initial state, @P(X0)@. If nothing then considered @(0.5,0.5)@
67
+
-> TransitionModel s R -- ^ Transition Model as @sxs@ matrix, @s@ being the number of states
68
+
-> V t b -- ^ evidence values vector map
69
+
-> V t (V s R) -- ^ a @txs@ matrix for each evidence value @t@, a vector capturing conditional probabilities for each state @s@
70
+
-> HMM s t R b
71
+
mkHMM_ mp xs evs = let p = fromMaybe (V $ Vec.replicate (dim xs) 0.5) mp
72
+
in
73
+
HMM p xs evs . (scaled <$>)
74
+
75
+
-- | A default HMM where both transition model states and evidence variables have boolean support
76
+
mkHMM :: KnownNat s => TransitionModel s R -> V 2 (V s R) -> HMM s 2 R Bool
77
+
mkHMM ts = mkHMM_ Nothing ts (V $ Vec.fromList [True, False])
78
+
79
+
-- | A default HMM with one state variable having boolean support
80
+
mkHMM1 :: V2 (V2 R) -> V2 (V2 R) -> HMM 2 2 R Bool
81
+
mkHMM1 ts ss = let f = (_V #) . over mapped toV
82
+
in
83
+
uncurry mkHMM $ over both f (ts , ss)
84
+
85
+
-- | create a @sxs@ diagonal matrix with corresponding posterior probabilities from the sensor model of the @HMM@ for an input sensor value
86
+
sensorDiagonal :: (KnownNat s, KnownNat t, Eq b)
87
+
=> HMM s t a b
88
+
-> b -- ^ sensor value
89
+
-> Maybe (SensorDiagonal s a)
90
+
sensorDiagonal hmm e = findIndexOf (sDist.folded) (== e) hmm >>= (\i -> hmm ^? sModel . ix i)
91
+
92
+
93
+
-- * Utility Functions
37
94
38
95
-- | Multiply each number by a constant such that the sum is 1.0
39
96
normalise :: (Fractional a, Foldable f, Functor f) => f a -> f a
···
44
101
in
45
102
(/ s) <$> xs
46
103
47
-
toHM :: (KnownNat s, KnownNat t) => V s (V t R) -> HM.L s t
48
-
toHM = HM.matrix . foldMap (Vec.toList . toVector)
49
-
50
-
fromHM :: (KnownNat s, KnownNat t) => HM.L s t -> V s (V t R)
51
-
fromHM m = V $ V <$> Vec.fromList (fmap Vec.fromList $ HMat.toLists $ HM.extract m)
52
-
104
+
-- | inverse a square matrix
53
105
inverse :: forall s. (KnownNat s) => V s (V s R) -> V s (V s R)
54
106
inverse = fromHM . HM.inv . toHM
55
107
108
+
-- | reverse contents of a 'Linear.V' vector
56
109
reverseV :: forall s a. (KnownNat s) => V s a -> V s a
57
110
reverseV = V . Vec.fromList . Papa.reverse . toList
58
111
59
112
extract1 :: (KnownNat s) => V s (V 1 a) -> V s a
60
113
extract1 = V . foldMap toVector
61
114
115
+
-- | Pointwise product of 2 @mxn@ 'V' vectors
62
116
infix 8 !**!
63
117
(!**!) :: (KnownNat m, KnownNat n) => V m (V n R) -> V m (V n R) -> V m (V n R)
64
118
as !**! bs = fromHM $ toHM as * toHM bs
65
119
120
+
-- | an unit vector having @s@ columns
66
121
unitColumn :: forall s. KnownNat s => Message s R
67
122
unitColumn = V $ Vec.replicate (fromIntegral $ natVal (Proxy :: Proxy s)) (toV $ V1 1.0)
68
123
124
+
-- ** Converting to/from __HMatrix__ matrices
69
125
70
-
-- | Hidden Markov Models
126
+
-- | take a @sxt@ 'Linear.V.V' matrix and give corresponding 'Numeric.LinearAlgebra.Static.L' version
127
+
toHM :: (KnownNat s, KnownNat t) => V s (V t R) -> HM.L s t
128
+
toHM = HM.matrix . foldMap (Vec.toList . toVector)
71
129
72
-
data HMM (s :: Nat) (t :: Nat) a b = HMM {
73
-
-- prior distribution which is used as initial forward message
74
-
_prior :: Message s a,
75
-
-- transition model with `s` states as `sxs` matrix
76
-
_tModel :: TransitionModel s a,
77
-
-- evidence value vector, each index maps to the corresponding evidence values
78
-
_sDist :: V t b,
79
-
-- sensor model with `t` evidence values having `sxs` diagonal matrix capturing their ditributions for each state `s`
80
-
_sModel :: SensorModel t s a
81
-
} deriving (Show)
82
-
makeLenses ''HMM
130
+
-- | take a @sxt@ 'Numeric.LinearAlgebra.Static.L' matrix and give corresponding 'Linear.V.V' version
131
+
fromHM :: (KnownNat s, KnownNat t) => HM.L s t -> V s (V t R)
132
+
fromHM m = V $ V <$> Vec.fromList (fmap Vec.fromList $ HMat.toLists $ HM.extract m)
83
133
84
134
85
-
-- | Smart HMM constructors
135
+
-- * Inference Algorithms
86
136
87
-
-- | Make HMM
88
-
-- `mp :: Maybe (Message s a)` -- Prior distribution on the initial state, `P(X0)`. If nothing then considered `(0.5,0.5)`
89
-
-- `xs :: TransitionModel s a` -- Transition Model as `sxs` matrix, `s` being the number of states
90
-
-- `evs :: V t b` - evidence values vector map
91
-
-- `es :: V t (V s a)` -- a `txs` matrix for each evidence value `t`, a vector capturing conditional probabilities for each state `s`
92
-
mkHMM_ :: (KnownNat s, KnownNat t) => Maybe (Message s R) -> TransitionModel s R -> V t b -> V t (V s R) -> HMM s t R b
93
-
mkHMM_ mp xs evs = let p = fromMaybe (V $ Vec.replicate (dim xs) 0.5) mp
94
-
in
95
-
HMM p xs evs . (scaled <$>)
96
-
97
-
-- | A default HMM where both transition model states and evidence variables have boolean support
98
-
mkHMM :: KnownNat s => TransitionModel s R -> V 2 (V s R) -> HMM s 2 R Bool
99
-
mkHMM ts = mkHMM_ Nothing ts (V $ Vec.fromList [True, False])
100
-
101
-
-- | A default HMM with one state variable having boolean support
102
-
mkHMM1 :: V2 (V2 R) -> V2 (V2 R) -> HMM 2 2 R Bool
103
-
mkHMM1 ts ss = let f = (_V #) . over mapped toV
104
-
in
105
-
uncurry mkHMM $ over both f (ts , ss)
106
-
107
-
sensorDiagonal :: (KnownNat s, KnownNat t, Eq b) => HMM s t a b -> b -> Maybe (SensorDiagonal s a)
108
-
sensorDiagonal hmm e = findIndexOf (sDist.folded) (== e) hmm >>= (\i -> hmm ^? sModel . ix i)
137
+
-- ** Forward-Backward
109
138
110
139
-- | Filtering message propagated forward
111
140
forward :: (KnownNat s) => HMM s t R b -> SensorDiagonal s R -> Message s R -> Message s R
···
115
144
backward :: (KnownNat s) => HMM s t R b -> SensorDiagonal s R -> Message s R -> Message s R
116
145
backward hmm ok1 bk2 = (hmm ^. tModel) !*! ok1 !*! bk2
117
146
147
+
-- | The forward–backward algorithm for smoothing: computing posterior prob- abilities of a sequence of states given a sequence of observations
148
+
forwardBackward :: forall s t u b . (KnownNat t, KnownNat s, KnownNat u, Eq b)
149
+
=> HMM s u R b -- ^ HMM model as a way to implement
150
+
-> V t b -- ^ list of evidences for each time step
151
+
-> V t (Distribution s R)
152
+
forwardBackward hmm evs = let
153
+
-- reifying the number of evidences
154
+
tNat = dim evs
155
+
-- forward messages from time t .. 0
156
+
fv :: V (t + 1) (Message s R)
157
+
fv = V $ Vec.fromList $ foldl' (\m@(x:_) e -> forward hmm (sensorDiagonal hmm e ^?! _Just) x : m) [hmm ^. prior] evs
158
+
-- getting rid of the prior message from the end of the list
159
+
-- so now forward messages vector is from time t .. 1
160
+
fv_0 :: V t (Message s R)
161
+
fv_0 = V $ Vec.fromList $ fv ^.. taking tNat traversed
118
162
119
-
-- | Fixed Lag Smoothing
163
+
-- backward messages from time 1 .. t
164
+
bs :: V t (Message s R)
165
+
bs = V $ Vec.fromList $ foldl' (\m@(x:_) e -> backward hmm (sensorDiagonal hmm e ^?! _Just) x : m) [unitColumn] $ reverseV evs
166
+
-- reversing the backward messages
167
+
-- so now the backward messages vector is from time t .. 1
168
+
revBs :: V t (Message s R)
169
+
revBs = reverseV bs
170
+
in
171
+
-- smoothing probabilities in reverse order of the list of evidences
172
+
-- i.e. starting from time t .. 1
173
+
liftA2 (\f b' -> extract1 $ normalise $ f !**! b') fv_0 revBs
174
+
120
175
176
+
-- ** Fixed Lag Smoothing
177
+
178
+
-- *** Data Types
121
179
-- | Persistent State
122
180
data Persistent (s :: Nat) a b d= Persistent {
123
-
_t :: d,
124
-
_f_msg :: Message s a,
125
-
_b :: V s (V s a),
126
-
_e_td_t :: Vec.Vector b
181
+
_t :: d, -- ^ current time
182
+
_f_msg :: Message s a, -- ^ the forward message @P(Xt|e1:t)@
183
+
_b :: V s (V s a), -- ^ the @d@-step backward transformation matrix
184
+
_e_td_t :: Vec.Vector b -- ^ double-ended list of evidence from @t − d@ to @t@
127
185
} deriving (Show)
186
+
187
+
-- |
188
+
-- === __Persistent lenses__
189
+
--
128
190
makeLenses ''Persistent
129
191
130
-
-- | State to start with
192
+
-- *** Smart constructor
193
+
-- | Initial persistent state where
194
+
--
195
+
-- * t = 1
196
+
--
197
+
-- * f_msg = hmm.prior
198
+
--
199
+
-- * b = identity matrix
200
+
--
201
+
-- * e_td_t = empty vector
131
202
persistentInit :: forall s b d. (KnownNat s, Integral d) => Message s R -> Persistent s R b d
132
203
persistentInit p = Persistent { _t = 1, _f_msg = p, _b = identity, _e_td_t = Vec.empty}
133
204
134
-
-- | Online Algorithm for smoothing with a fixed time lag of `d` steps
135
-
-- hmm -- HMM model
136
-
-- d -- length of lag
137
-
-- e -- evidence at time t
138
-
fixedLagSmoothing :: forall s u b d. (KnownNat s, KnownNat u, Eq b, Integral d) => HMM s u R b -> d -> b -> MaybeT (State (Persistent s R b d)) (Distribution s R)
205
+
-- *** Algo
206
+
-- | Online Algorithm for smoothing with a fixed time lag of @d@ steps
207
+
fixedLagSmoothing :: forall s u b d. (KnownNat s, KnownNat u, Eq b, Integral d)
208
+
=> HMM s u R b -- ^ HMM model
209
+
-> d -- ^ length of lag
210
+
-> b -- ^ evidence at time @t@
211
+
-> MaybeT (State (Persistent s R b d)) (Distribution s R)
139
212
fixedLagSmoothing hmm d e = do
140
213
e_td_t %= flip Vec.snoc e
141
214
···
159
232
t += 1
160
233
mzero
161
234
235
+
-- ** [/WIP/] Cost based HMM(cHMM)
236
+
-- Inference algo from cost-HMM paper
162
237
163
-
-- | The forward–backward algorithm for smoothing: computing posterior prob- abilities of a sequence of states given a sequence of observations
164
-
-- hmm - HMM model as a way to implement
165
-
-- evs - list of evidences for each time step
166
-
forwardBackward :: forall s t u b . (KnownNat t, KnownNat s, KnownNat u, Eq b) => HMM s u R b -> V t b -> V t (Distribution s R)
167
-
forwardBackward hmm evs = let
168
-
-- reifying the number of evidences
169
-
tNat = dim evs
170
-
-- forward messages from time t .. 0
171
-
fv :: V (t + 1) (Message s R)
172
-
fv = V $ Vec.fromList $ foldl' (\m@(x:_) e -> forward hmm (sensorDiagonal hmm e ^?! _Just) x : m) [hmm ^. prior] evs
173
-
-- getting rid of the prior message from the end of the list
174
-
-- so now forward messages vector is from time t .. 1
175
-
fv_0 :: V t (Message s R)
176
-
fv_0 = V $ Vec.fromList $ fv ^.. taking tNat traversed
238
+
cHMMInfer :: forall s u t b i. (KnownNat s, KnownNat u, KnownNat t, Eq b, Ord b, Integral i, Show i, Show b)
239
+
=> HMM s u R b -- ^ HMM model
240
+
-> V t i -- ^ state trajectory for time @1 .. t@; represented as @0@-based index of state
241
+
-> V t b -- ^ evidence/observation trajectory for time @1 .. t@
242
+
-> IO R
243
+
cHMMInfer hmm xs evs = do
244
+
let r = (\x y z -> (x * y) / z) <$> y0 <*> y1 <*> y2 $ zs
245
+
print zs
246
+
print $ y0 zs
247
+
print $ y1 zs
248
+
return r
249
+
where
250
+
y0 :: Vec.Vector (i, b) -> R
251
+
y0 = foldl' (\m (s,o') ->
252
+
m *
253
+
sensorDiagonal hmm o' ^?! _Just . ix (fromIntegral s) . ix (fromIntegral s))
254
+
1
177
255
178
-
-- backward messages from time 1 .. t
179
-
bs :: V t (Message s R)
180
-
bs = V $ Vec.fromList $ foldl' (\m@(x:_) e -> backward hmm (sensorDiagonal hmm e ^?! _Just) x : m) [unitColumn] $ reverseV evs
181
-
-- reversing the backward messages
182
-
-- so now the backward messages vector is from time t .. 1
183
-
revBs :: V t (Message s R)
184
-
revBs = reverseV bs
185
-
in
186
-
-- smoothing probabilities in reverse order of the list of evidences
187
-
-- i.e. starting from time t .. 1
188
-
liftA2 (\f b' -> extract1 $ normalise $ f !**! b') fv_0 revBs
256
+
y1 :: Vec.Vector (i, b) -> R
257
+
y1 vs = ifoldl (\idx m (s,_) ->
258
+
m *
259
+
bool
260
+
(hmm ^?! tModel . ix (vs ^?! ix (idx - 1) . _1 . to fromIntegral) . ix (fromIntegral s))
261
+
(((Linear.transpose $ hmm ^. prior) !*! hmm ^. tModel) ^?! ix 0 . ix (fromIntegral s))
262
+
(idx == 0))
263
+
1 vs
264
+
265
+
y2 :: KnownNat s => Vec.Vector (i, b) -> R
266
+
y2 vs = sum . extract1 $
267
+
ifoldl (\idx m (_,o') ->
268
+
bool
269
+
((sensorDiagonal hmm o' ^?! _Just) !*! Linear.transpose (hmm ^. tModel) !*! m)
270
+
((sensorDiagonal hmm o' ^?! _Just) !*! Linear.transpose (hmm ^. tModel) !*! hmm ^. prior)
271
+
(idx == 0)
272
+
)
273
+
(unitColumn :: Message s R)
274
+
vs
275
+
276
+
zs :: Vec.Vector (i, b)
277
+
zs = Vec.zip (toVector xs) (toVector evs)
278
+
279
+
-- * Sample HMM Models
189
280
190
-
-- | execution
191
281
umbrellaHMM :: HMM 2 2 R Bool
192
282
umbrellaHMM = mkHMM1 (V2 (V2 0.7 0.3) (V2 0.3 0.7)) (V2 (V2 0.9 0.2) (V2 0.1 0.8))
283
+
284
+
swarmRobotsHMM :: HMM 4 4 R Int
285
+
swarmRobotsHMM = mkHMM_
286
+
(Just $ V $ Vec.replicate 4 (toV $ V1 0.25))
287
+
(V $ Vec.fromList
288
+
[
289
+
(V $ Vec.fromList [0.7393859 , 0.2233470 , 0.0330670 , 0.0042001]),
290
+
(V $ Vec.fromList [0.0928202 , 0.7172515 , 0.1719556 , 0.0179727]),
291
+
(V $ Vec.fromList [0.0139341 , 0.1755173 , 0.6936519 , 0.1168966]),
292
+
(V $ Vec.fromList [0.0047083 , 0.0403520 , 0.2561893 , 0.6987504])
293
+
])
294
+
(V $ Vec.fromList [1..4])
295
+
(Linear.transpose $ (V $ Vec.fromList
296
+
[
297
+
(V $ Vec.fromList [1.00000, 0.00000, 0.00000, 0.00000]),
298
+
(V $ Vec.fromList [0.10000, 0.90000, 0.00000, 0.00000]),
299
+
(V $ Vec.fromList [0.01000, 0.18000, 0.81000, 0.00000]),
300
+
(V $ Vec.fromList [0.00100, 0.02700, 0.24300, 0.72900])
301
+
]))
193
302
303
+
-- * Main
194
304
runFLSAlgo :: [Bool] -> Integer -> [Maybe (Distribution 2 String)]
195
305
runFLSAlgo bs d = let hmm = mkHMM1 (V2 (V2 0.7 0.3) (V2 0.3 0.7)) (V2 (V2 0.9 0.2) (V2 0.1 0.8))
196
306
initState = persistentInit (V $ Vec.fromList [V $ Vec.singleton 0.5, V $ Vec.singleton 0.5]) :: Persistent 2 R Bool Integer
···
202
312
203
313
204
314
runFBAlgo :: [Bool] -> [Distribution 2 R]
205
-
runFBAlgo bs = reifyVectorNat (Vec.fromList bs) (toList . forwardBackward umbrellaHMM)
315
+
runFBAlgo bs = reifyVectorNat (Vec.fromList bs) (toList . forwardBackward umbrellaHMM)