Haskell implementations of Hidden Markov Models & related algorithms from AIMA book
documentations
+121
-84
+1
-1
README.md
+1
-1
README.md
···
1
-
Haskell implementation of Hidden Markov Models & related algorithms from AIMA book
1
+
Haskell implementations of Hidden Markov Models & related algorithms from AIMA book and using them for approaches mentioned in the paper [Assessing the resilience of stochastic dynamic systems under partial observability](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0202337)
+93
-1
src/cHMM.hs
+93
-1
src/cHMM.hs
···
10
10
-- Maintainer : Kaushik Chakraborty <git@kaushikc.org>
11
11
-- Stability : experimental
12
12
--
13
+
-- Cost based HMM
14
+
-- https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0202337
13
15
---------------------------------
14
16
15
17
module CHMM where
···
18
20
import GHC.TypeNats
19
21
import Data.Proxy (Proxy(..))
20
22
import Linear.V
23
+
import Linear
21
24
import Combinatorics
22
25
import qualified Data.Vector as Vec
23
26
24
27
import HMM
25
28
-----------------------
26
29
27
-
trajectory :: forall n t a . (KnownNat n, KnownNat t) => Proxy t -> V n a -> V (n ^ t) (V n a)
30
+
-- * Documentation
31
+
32
+
-- | For a specific resiliency property that need evaluation, following algorithm infer the probability of the same property in linear time taking into account corresponding observations and state trajectories
33
+
cHMMInfer :: forall s u t b i. (KnownNat s, KnownNat u, KnownNat t, Eq b, Ord b, Integral i, Show i, Show b)
34
+
=> HMM s u R b -- ^ HMM model
35
+
-> V t b -- ^ evidence/observation trajectory for time @1 .. t@
36
+
-> V t i -- ^ state trajectory for time @1 .. t@; represented as @0@-based index of state
37
+
-> R -- ^ property probability
38
+
cHMMInfer hmm evs xs = (\x y z -> (x * y) / z) <$> y0 <*> y1 <*> y2 $ zs
39
+
where
40
+
y0 :: Vec.Vector (i, b) -> R
41
+
y0 = foldl' (\m (s,o') ->
42
+
m *
43
+
sensorDiagonal hmm o' ^?! _Just . ix (fromIntegral s) . ix (fromIntegral s))
44
+
1
45
+
46
+
y1 :: Vec.Vector (i, b) -> R
47
+
y1 vs = ifoldl (\idx m (s,_) ->
48
+
m *
49
+
bool
50
+
(hmm ^?! tModel . ix (vs ^?! ix (idx - 1) . _1 . to fromIntegral) . ix (fromIntegral s))
51
+
((Linear.transpose (hmm ^. prior) !*! hmm ^. tModel) ^?! ix 0 . ix (fromIntegral s))
52
+
(idx == 0))
53
+
1 vs
54
+
55
+
y2 :: KnownNat s => Vec.Vector (i, b) -> R
56
+
y2 vs = sum . extract1 $
57
+
ifoldl (\idx m (_,o') ->
58
+
bool
59
+
((sensorDiagonal hmm o' ^?! _Just) !*! Linear.transpose (hmm ^. tModel) !*! m)
60
+
((sensorDiagonal hmm o' ^?! _Just) !*! Linear.transpose (hmm ^. tModel) !*! hmm ^. prior)
61
+
(idx == 0)
62
+
)
63
+
(unitColumn :: Message s R)
64
+
vs
65
+
66
+
zs :: Vec.Vector (i, b)
67
+
zs = Vec.zip (toVector xs) (toVector evs)
68
+
69
+
-- | a trajectory is a sequence of variables which can be random variables (r.v) as part of either sensor or transition model of a HMM or values from the support of some r.v.
70
+
trajectory :: forall n t a . (KnownNat n, KnownNat t)
71
+
=> Proxy t -- ^ choice parameter
72
+
-> V n a -- ^ superset of choice variables which is divided into sets each of length provided by the choice parameter
73
+
-> V (n ^ t) (V t a) -- ^ final set with @n^t@ ordered subsets each of length @t@
28
74
trajectory p = V . Vec.fromList . fmap (V . Vec.fromList) . variateRep t' . toList
29
75
where
30
76
t' = fromIntegral $ natVal p
31
77
78
+
-- | HMM model representing the swarm robotics use case of the paper
79
+
--
80
+
-- Transition Model
81
+
--
82
+
-- @
83
+
-- [0.7393859 , 0.2233470 , 0.0330670 , 0.0042001]
84
+
-- [0.0928202 , 0.7172515 , 0.1719556 , 0.0179727]
85
+
-- [0.0139341 , 0.1755173 , 0.6936519 , 0.1168966]
86
+
-- [0.0047083 , 0.0403520 , 0.2561893 , 0.6987504]
87
+
-- @
88
+
--
89
+
-- Sensor Model
90
+
--
91
+
-- @
92
+
-- [1.00000, 0.00000, 0.00000, 0.00000]
93
+
-- [0.10000, 0.90000, 0.00000, 0.00000]
94
+
-- [0.01000, 0.18000, 0.81000, 0.00000]
95
+
-- [0.00100, 0.02700, 0.24300, 0.72900]
96
+
-- @
97
+
--
98
+
swarmRobotsHMM :: HMM 4 4 R Int
99
+
swarmRobotsHMM = mkHMM_
100
+
(Just $ V $ Vec.replicate 4 (toV $ V1 0.25))
101
+
(V $ Vec.fromList
102
+
[
103
+
V $ Vec.fromList [0.7393859 , 0.2233470 , 0.0330670 , 0.0042001],
104
+
V $ Vec.fromList [0.0928202 , 0.7172515 , 0.1719556 , 0.0179727],
105
+
V $ Vec.fromList [0.0139341 , 0.1755173 , 0.6936519 , 0.1168966],
106
+
V $ Vec.fromList [0.0047083 , 0.0403520 , 0.2561893 , 0.6987504]
107
+
])
108
+
(V $ Vec.fromList [1..4])
109
+
(Linear.transpose (V $ Vec.fromList
110
+
[
111
+
V $ Vec.fromList [1.00000, 0.00000, 0.00000, 0.00000],
112
+
V $ Vec.fromList [0.10000, 0.90000, 0.00000, 0.00000],
113
+
V $ Vec.fromList [0.01000, 0.18000, 0.81000, 0.00000],
114
+
V $ Vec.fromList [0.00100, 0.02700, 0.24300, 0.72900]
115
+
]))
116
+
117
+
-- | inferring the /l-resistance/ property of the swarm robots as represented in the 'swarmRobotsHMM' model given a series of observation and a paticular /l/ value to check against
118
+
--
119
+
-- E.g.
120
+
--
121
+
-- >>> runCHMMInfer [3,2,3] 2
122
+
-- 0.7756162284376005
123
+
--
32
124
runCHMMInfer :: [Int] -> Int -> R
33
125
runCHMMInfer evs l = let v = Vec.fromList [0 .. dim (swarmRobotsHMM ^. tModel) - 1]
34
126
n = length evs
+27
-82
src/hmm.hs
+27
-82
src/hmm.hs
···
13
13
--
14
14
---------------------------------
15
15
16
-
module HMM where
16
+
module HMM
17
+
(
18
+
HMM(..),
19
+
-- ** Smart Constructors
20
+
mkHMM_, mkHMM, mkHMM1, sensorDiagonal,
21
+
-- * Utility Functions
22
+
normalise, inverse, reverseV, extract1, (!**!), unitColumn,
23
+
-- ** Converting to/from HMatrix matrices
24
+
toHM, fromHM,
25
+
-- * Inference Algorithms
26
+
-- ** Forward-Backward
27
+
forward, backward, forwardBackward,
28
+
-- ** Fixed Lag Smoothing
29
+
fixedLagSmoothing,
30
+
-- * Sample HMM Models
31
+
umbrellaHMM,
32
+
-- * Miscellaneous
33
+
-- ** Type Synonyms
34
+
R, TransitionModel, SensorDiagonal, SensorModel, Message, Distribution,
35
+
-- ** Data Types
36
+
Persistent(..),
37
+
-- ** Test Functions
38
+
runFLSAlgo, runFBAlgo,
39
+
-- ** Lenses
40
+
prior, tModel, sDist, sModel
41
+
) where
17
42
18
43
import Text.Printf (printf)
19
44
···
33
58
import qualified Numeric.LinearAlgebra as HMat
34
59
-----------------------------------------------------------------------------------
35
60
36
-
-- * Type synonyms
37
61
type R = Double
38
62
type TransitionModel (s :: Nat) a = V s (V s a)
39
63
type SensorDiagonal (s :: Nat) a = V s (V s a)
···
41
65
type Message (s :: Nat) a = V s (V 1 a)
42
66
type Distribution (s :: Nat) a = V s a
43
67
44
-
-- * HMM
45
-
46
68
-- | Hidden Markov Models
47
69
data HMM (s :: Nat) (t :: Nat) a b = HMM {
48
70
-- | prior distribution which is used as initial forward message
···
55
77
_sModel :: SensorModel t s a
56
78
} deriving (Show)
57
79
58
-
-- |
59
-
-- === __HMM lenses__
60
-
--
61
80
makeLenses ''HMM
62
81
63
-
-- ** Smart constructors
64
82
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)@
65
83
-> TransitionModel s R -- ^ Transition Model as @sxs@ matrix, @s@ being the number of states
66
84
-> V t b -- ^ evidence values vector map
···
88
106
sensorDiagonal hmm e = findIndexOf (sDist.folded) (== e) hmm >>= (\i -> hmm ^? sModel . ix i)
89
107
90
108
91
-
-- * Utility Functions
92
109
93
110
-- | Multiply each number by a constant such that the sum is 1.0
94
111
normalise :: (Fractional a, Foldable f, Functor f) => f a -> f a
···
119
136
unitColumn :: forall s. KnownNat s => Message s R
120
137
unitColumn = V $ Vec.replicate (fromIntegral $ natVal (Proxy :: Proxy s)) (toV $ V1 1.0)
121
138
122
-
-- ** Converting to/from __HMatrix__ matrices
123
-
124
139
-- | take a @sxt@ 'Linear.V.V' matrix and give corresponding 'Numeric.LinearAlgebra.Static.L' version
125
140
toHM :: (KnownNat s, KnownNat t) => V s (V t R) -> HM.L s t
126
141
toHM = HM.matrix . foldMap (Vec.toList . toVector)
···
128
143
-- | take a @sxt@ 'Numeric.LinearAlgebra.Static.L' matrix and give corresponding 'Linear.V.V' version
129
144
fromHM :: (KnownNat s, KnownNat t) => HM.L s t -> V s (V t R)
130
145
fromHM m = V $ V <$> Vec.fromList (fmap Vec.fromList $ HMat.toLists $ HM.extract m)
131
-
132
-
133
-
-- * Inference Algorithms
134
-
135
-
-- ** Forward-Backward
136
146
137
147
-- | Filtering message propagated forward
138
148
forward :: (KnownNat s) => HMM s t R b -> SensorDiagonal s R -> Message s R -> Message s R
···
171
181
liftA2 (\f b' -> extract1 $ normalise $ f !**! b') fv_0 revBs
172
182
173
183
174
-
-- ** Fixed Lag Smoothing
175
-
176
-
-- *** Data Types
177
184
-- | Persistent State
178
185
data Persistent (s :: Nat) a b d= Persistent {
179
186
_t :: d, -- ^ current time
···
182
189
_e_td_t :: Vec.Vector b -- ^ double-ended list of evidence from @t − d@ to @t@
183
190
} deriving (Show)
184
191
185
-
-- |
186
-
-- === __Persistent lenses__
187
-
--
188
192
makeLenses ''Persistent
189
193
190
-
-- *** Smart constructor
191
194
-- | Initial persistent state where
192
195
--
193
196
-- * t = 1
···
200
203
persistentInit :: forall s b d. (KnownNat s, Integral d) => Message s R -> Persistent s R b d
201
204
persistentInit p = Persistent { _t = 1, _f_msg = p, _b = identity, _e_td_t = Vec.empty}
202
205
203
-
-- *** Algo
204
206
-- | Online Algorithm for smoothing with a fixed time lag of @d@ steps
205
207
fixedLagSmoothing :: forall s u b d. (KnownNat s, KnownNat u, Eq b, Integral d)
206
208
=> HMM s u R b -- ^ HMM model
···
230
232
t += 1
231
233
mzero
232
234
233
-
-- ** Cost based HMM(cHMM)
234
-
-- Inference algo from cost-HMM paper
235
-
cHMMInfer :: forall s u t b i. (KnownNat s, KnownNat u, KnownNat t, Eq b, Ord b, Integral i, Show i, Show b)
236
-
=> HMM s u R b -- ^ HMM model
237
-
-> V t b -- ^ evidence/observation trajectory for time @1 .. t@
238
-
-> V t i -- ^ state trajectory for time @1 .. t@; represented as @0@-based index of state
239
-
-> R
240
-
cHMMInfer hmm evs xs = (\x y z -> (x * y) / z) <$> y0 <*> y1 <*> y2 $ zs
241
-
where
242
-
y0 :: Vec.Vector (i, b) -> R
243
-
y0 = foldl' (\m (s,o') ->
244
-
m *
245
-
sensorDiagonal hmm o' ^?! _Just . ix (fromIntegral s) . ix (fromIntegral s))
246
-
1
247
235
248
-
y1 :: Vec.Vector (i, b) -> R
249
-
y1 vs = ifoldl (\idx m (s,_) ->
250
-
m *
251
-
bool
252
-
(hmm ^?! tModel . ix (vs ^?! ix (idx - 1) . _1 . to fromIntegral) . ix (fromIntegral s))
253
-
((Linear.transpose (hmm ^. prior) !*! hmm ^. tModel) ^?! ix 0 . ix (fromIntegral s))
254
-
(idx == 0))
255
-
1 vs
256
-
257
-
y2 :: KnownNat s => Vec.Vector (i, b) -> R
258
-
y2 vs = sum . extract1 $
259
-
ifoldl (\idx m (_,o') ->
260
-
bool
261
-
((sensorDiagonal hmm o' ^?! _Just) !*! Linear.transpose (hmm ^. tModel) !*! m)
262
-
((sensorDiagonal hmm o' ^?! _Just) !*! Linear.transpose (hmm ^. tModel) !*! hmm ^. prior)
263
-
(idx == 0)
264
-
)
265
-
(unitColumn :: Message s R)
266
-
vs
267
-
268
-
zs :: Vec.Vector (i, b)
269
-
zs = Vec.zip (toVector xs) (toVector evs)
270
-
271
-
-- * Sample HMM Models
272
236
umbrellaHMM :: HMM 2 2 R Bool
273
237
umbrellaHMM = mkHMM1 (V2 (V2 0.7 0.3) (V2 0.3 0.7)) (V2 (V2 0.9 0.2) (V2 0.1 0.8))
274
238
275
-
swarmRobotsHMM :: HMM 4 4 R Int
276
-
swarmRobotsHMM = mkHMM_
277
-
(Just $ V $ Vec.replicate 4 (toV $ V1 0.25))
278
-
(V $ Vec.fromList
279
-
[
280
-
V $ Vec.fromList [0.7393859 , 0.2233470 , 0.0330670 , 0.0042001],
281
-
V $ Vec.fromList [0.0928202 , 0.7172515 , 0.1719556 , 0.0179727],
282
-
V $ Vec.fromList [0.0139341 , 0.1755173 , 0.6936519 , 0.1168966],
283
-
V $ Vec.fromList [0.0047083 , 0.0403520 , 0.2561893 , 0.6987504]
284
-
])
285
-
(V $ Vec.fromList [1..4])
286
-
(Linear.transpose (V $ Vec.fromList
287
-
[
288
-
V $ Vec.fromList [1.00000, 0.00000, 0.00000, 0.00000],
289
-
V $ Vec.fromList [0.10000, 0.90000, 0.00000, 0.00000],
290
-
V $ Vec.fromList [0.01000, 0.18000, 0.81000, 0.00000],
291
-
V $ Vec.fromList [0.00100, 0.02700, 0.24300, 0.72900]
292
-
]))
293
239
294
-
-- * Main
295
240
runFLSAlgo :: [Bool] -> Integer -> [Maybe (Distribution 2 String)]
296
241
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))
297
242
initState = persistentInit (V $ Vec.fromList [V $ Vec.singleton 0.5, V $ Vec.singleton 0.5]) :: Persistent 2 R Bool Integer
···
303
248
304
249
305
250
runFBAlgo :: [Bool] -> [Distribution 2 R]
306
-
runFBAlgo bs = reifyVectorNat (Vec.fromList bs) (toList . forwardBackward umbrellaHMM)
251
+
runFBAlgo bs = reifyVectorNat (Vec.fromList bs) (toList . forwardBackward umbrellaHMM)