packages feed

learning-hmm 0.1.1.1 → 0.2.0.0

raw patch · 12 files changed

+382/−365 lines, 12 filesdep +deepseqdep −logfloatPVP ok

version bump matches the API change (PVP)

Dependencies added: deepseq

Dependencies removed: logfloat

API changes (from Hackage documentation)

- Learning.HMM: init :: (Ord s, Ord o) => [s] -> [o] -> RVar (HMM s o)
+ Learning.HMM: init :: (Eq s, Eq o) => [s] -> [o] -> RVar (HMM s o)
- Learning.HMM: withEmission :: (Ord s, Ord o) => HMM s o -> [o] -> HMM s o
+ Learning.HMM: withEmission :: (Eq s, Eq o) => HMM s o -> [o] -> HMM s o

Files

CHANGES.md view
@@ -1,6 +1,10 @@ Revision history for Haskell package learning-hmm === +## Version 0.2.0.0+- Remove dependency on the 'logfloat' package+- Performance improvements+ ## Version 0.1.1.0 - Add function `init` for random initialization - Add function `simulate` for running a Markov process
learning-hmm.cabal view
@@ -1,5 +1,5 @@ name:                learning-hmm-version:             0.1.1.1+version:             0.2.0.0 stability:           experimental  synopsis:            Yet another library for hidden Markov models@@ -27,14 +27,13 @@   exposed-modules:   Learning.HMM   other-modules:     Data.Random.Distribution.Categorical.Util                    , Data.Random.Distribution.Simplex-                   , Data.Random.Distribution.Uniform.Util-                   , Data.Vector.Util-                   , Data.Vector.Util.LinearAlgebra+                   , Data.Vector.Generic.Util+                   , Data.Vector.Generic.Util.LinearAlgebra                    , Learning.HMM.Internal   -- other-extensions:     build-depends:     base >=4.7 && <4.8                    , containers-                   , logfloat+                   , deepseq                    , random-fu                    , random-source                    , vector
src/Data/Random/Distribution/Categorical/Util.hs view
@@ -2,11 +2,10 @@  module Data.Random.Distribution.Categorical.Util () where +import Data.Maybe (fromMaybe) import Data.Random.Distribution (PDF, pdf)-import Data.Random.Distribution.Categorical (Categorical, toList, totalWeight)+import Data.Random.Distribution.Categorical (Categorical, toList) import Data.Tuple (swap)  instance Eq a => PDF (Categorical Double) a where-  pdf cat a = case lookup a $ map swap $ toList cat of-                Nothing -> 0-                Just p  -> p / totalWeight cat+  pdf cat a = fromMaybe 0 (lookup a $ map swap $ toList cat)
src/Data/Random/Distribution/Simplex.hs view
@@ -11,7 +11,6 @@     , fractionalStdSimplex     ) where -import Control.Applicative import Control.Monad import Data.List import Data.Random.RVar
− src/Data/Random/Distribution/Uniform/Util.hs 
@@ -1,17 +0,0 @@-{-# LANGUAGE FlexibleContexts, FlexibleInstances, MultiParamTypeClasses #-}--module Data.Random.Distribution.Uniform.Util () where--import Control.Applicative ((<$>))-import Data.Number.LogFloat (LogFloat, logFloat, fromLogFloat)-import Data.Random.Distribution (Distribution)-import Data.Random.Distribution.Uniform -- (StdUniform(..), Uniform(..), doubleUniform)-import Data.Random (rvarT)-import Data.Random.Source (getRandomDouble)--instance Distribution Uniform LogFloat where-  rvarT (Uniform a b) = do x <- doubleUniform (fromLogFloat a) (fromLogFloat b)-                           return $ logFloat x--instance Distribution StdUniform LogFloat where-  rvarT _ = logFloat <$> getRandomDouble
+ src/Data/Vector/Generic/Util.hs view
@@ -0,0 +1,19 @@+-- | Miscellaneous utility functions for "Data.Vector"+module Data.Vector.Generic.Util (+    frequencies+  ) where++import Data.Map.Strict (Map)+import qualified Data.Map.Strict as M (empty, insertWith)+import Data.Vector.Generic (Vector, foldl')++-- $setup+-- >>> :module + Data.Vector++-- | @frequencies xs@ returns a 'Map' from distinct items in @xs@ to+-- the number of times they appear.+--+-- >>> frequencies $ fromList "bra bra bar"+-- fromList [(' ',2),('a',3),('b',3),('r',3)]+frequencies :: (Ord a, Vector v a, Num n) => v a -> Map a n+frequencies = foldl' (\m k -> M.insertWith (+) k 1 m) M.empty
+ src/Data/Vector/Generic/Util/LinearAlgebra.hs view
@@ -0,0 +1,141 @@+-- | Operators commonly used in the basic linear algebra. Note that all the+--   functions defined here do not check the dimension/length of+--   vectors/matrices.+module Data.Vector.Generic.Util.LinearAlgebra (+  -- * Pairwise operators+    (>+>)+  -- , (>->)+  , (>.>)+  , (>/>)+  , (#+#)+  -- , (#-#)+  -- , (#.#)+  -- , (#/#)++  -- * Scalar-vector/vector-scalar operators+  -- , (+>)+  -- , (->)+  , (.>)+  -- , (/>)+  -- , (>+)+  -- , (>-)+  -- , (>.)+  , (>/)++  -- * Scalar-matrix/matrix-scalar operators+  -- , (+#)+  -- , (-#)+  -- , (.#)+  -- , (/#)+  -- , (#+)+  -- , (#-)+  -- , (#.)+  , (#/)++  -- * Dot and matrix-vector/vector-matrix products+  , (<.>)+  , (#.>)+  , (<.#)++  -- * Unary operators+  , transpose+  ) where++import Prelude hiding (any, head, map, null, sum, tail, zipWith)+import Data.Vector.Generic (+    Vector, any, cons, convert, empty, head, map, null, sum, tail, zipWith+  )++-- $setup+-- >>> :module + Data.Vector++-- | Pairwise addition between two vectors+--+-- >>> fromList [1, 2] >+> fromList [3, 4 :: Int]+-- fromList [4,6]+(>+>) :: (Num a, Vector v a) => v a -> v a -> v a+{-# INLINE (>+>) #-}+u >+> v = zipWith (+) u v++-- | Pairwise product between two vectors+--+-- >>> fromList [1, 2] >.> fromList [3, 4 :: Double]+-- fromList [3.0,8.0]+(>.>) :: (Num a, Vector v a) => v a -> v a -> v a+{-# INLINE (>.>) #-}+u >.> v = zipWith (*) u v++-- | Pairwise division between two vectors+--+-- >>> fromList [1, 2] >/> fromList [3, 4 :: Double]+-- fromList [0.3333333333333333,0.5]+(>/>) :: (Fractional a, Vector v a) => v a -> v a -> v a+{-# INLINE (>/>) #-}+u >/> v = zipWith (/) u v++-- | Pairwise addition between two matrices+--+-- >>> fromList [fromList [1, 2], fromList [3, 4]] #+# fromList [fromList [5, 6], fromList [7, 8 :: Int]]+-- fromList [fromList [6,8],fromList [10,12]]+(#+#) :: (Num a, Vector v a, Vector w (v a)) => w (v a) -> w (v a) -> w (v a)+{-# INLINE (#+#) #-}+m #+# n = zipWith (>+>) m n++-- | Scalar-vector product+--+-- >>> 2 .> fromList [1, 2 :: Integer]+-- fromList [2,4]+(.>) :: (Num a, Vector v a) => a -> v a -> v a+{-# INLINE (.>) #-}+s .> v = map (s *) v++-- | Vector-scalar division+--+-- >>> fromList [1, 2 :: Double] >/ 2+-- fromList [0.5,1.0]+(>/) :: (Fractional a, Vector v a) => v a -> a -> v a+{-# INLINE (>/) #-}+v >/ s = map (/ s) v++-- | Matrix-scalar division+--+-- >>> fromList [fromList [1, 2], fromList [3, 4 :: Double]] #/ 2+-- fromList [fromList [0.5,1.0],fromList [1.5,2.0]]+(#/) :: (Fractional a, Vector v a, Vector w (v a)) => w (v a) -> a -> w (v a)+{-# INLINE (#/) #-}+m #/ s = map (>/ s) m++-- | Dot product+--+-- >>> fromList [1, 2] <.> fromList [3, 4 :: Int]+-- 11+(<.>) :: (Num a, Vector v a) => v a -> v a -> a+{-# INLINE (<.>) #-}+u <.> v = sum $ u >.> v++-- | Matrix-vector product+--+-- >>> fromList [fromList [1, 2], fromList [3, 4]] #.> fromList [1, 2 :: Double]+-- fromList [5.0,11.0]+(#.>) :: (Num a, Vector v a, Vector w (v a), Vector w a) => w (v a) -> v a -> v a+{-# INLINE (#.>) #-}+m #.> v = convert $ map (<.> v) m++-- | Vector-matrix product+--+-- >>> fromList [1, 2 :: Double] <.# fromList [fromList [1, 2], fromList [3, 4]]+-- fromList [7.0,10.0]+(<.#) :: (Num a, Vector v a, Vector w (v a), Vector w a) => v a -> w (v a) -> v a+{-# INLINE (<.#) #-}+v <.# m | any null m = empty+        | otherwise  = (v <.> convert (map head m)) `cons` (v <.# map tail m)++-- | Matrix transpose+--+-- >>> transpose $ fromList [fromList "ab", fromList "cd"]+-- fromList [fromList "ac",fromList "bd"]+transpose :: (Vector v a, Vector w (v a), Vector w a) => w (v a) -> w (v a)+{-# INLINE transpose #-}+transpose m+  | any null m = empty+  | otherwise  = convert (map head m) `cons` transpose (map tail m)
− src/Data/Vector/Util.hs 
@@ -1,13 +0,0 @@--- | Miscellaneous utility functions for "Data.Vector"-module Data.Vector.Util (-    unsafeElemIndex-  ) where--import Data.Maybe (fromJust)-import Data.Vector (Vector, elemIndex)---- | Return the index of the first occurrence of the given element or throw---   an error if no such occurrence.-{-# INLINE unsafeElemIndex #-}-unsafeElemIndex :: Eq a => a -> Vector a -> Int-unsafeElemIndex e = fromJust . elemIndex e
− src/Data/Vector/Util/LinearAlgebra.hs 
@@ -1,139 +0,0 @@--- | Operators commonly used in the basic linear algebra. Note that all the---   functions defined here do not check the dimension/length of---   vectors/matrices.-module Data.Vector.Util.LinearAlgebra (-  -- * Pairwise operators-    (>+>)-  -- , (>->)-  , (>.>)-  , (>/>)-  , (#+#)-  -- , (#-#)-  -- , (#.#)-  -- , (#/#)--  -- * Scalar-vector/vector-scalar operators-  -- , (+>)-  -- , (->)-  , (.>)-  -- , (/>)-  -- , (>+)-  -- , (>-)-  -- , (>.)-  , (>/)--  -- * Scalar-matrix/matrix-scalar operators-  -- , (+#)-  -- , (-#)-  -- , (.#)-  -- , (/#)-  -- , (#+)-  -- , (#-)-  -- , (#.)-  , (#/)--  -- * Dot and matrix-vector/vector-matrix products-  , (<.>)-  , (#.>)-  , (<.#)--  -- * Unary operators-  , transpose-  ) where--import Prelude hiding (any, head, map, null, sum, tail, zipWith)-import Data.Vector (Vector, any, cons, empty, head, map, null, sum, tail, zipWith)---- $setup--- >>> :module + Data.Vector---- | Pairwise addition between two vectors------ >>> fromList [1, 2] >+> fromList [3, 4 :: Int]--- fromList [4,6]-(>+>) :: Num a => Vector a -> Vector a -> Vector a-{-# INLINE (>+>) #-}-u >+> v = zipWith (+) u v---- | Pairwise product between two vectors------ >>> fromList [1, 2] >.> fromList [3, 4 :: Double]--- fromList [3.0,8.0]-(>.>) :: Num a => Vector a -> Vector a -> Vector a-{-# INLINE (>.>) #-}-u >.> v = zipWith (*) u v---- | Pairwise division between two vectors------ >>> fromList [1, 2] >/> fromList [3, 4 :: Double]--- fromList [0.3333333333333333,0.5]-(>/>) :: Fractional a => Vector a -> Vector a -> Vector a-{-# INLINE (>/>) #-}-u >/> v = zipWith (/) u v---- | Pairwise addition between two matrices------ >>> fromList [fromList [1, 2], fromList [3, 4]] #+# fromList [fromList [5, 6], fromList [7, 8 :: Int]]--- fromList [fromList [6,8],fromList [10,12]]-(#+#) :: Num a => Vector (Vector a) -> Vector (Vector a) -> Vector (Vector a)-{-# INLINE (#+#) #-}-m #+# n = zipWith (>+>) m n---- | Scalar-vector product------ >>> 2 .> fromList [1, 2 :: Integer]--- fromList [2,4]-(.>) :: Num a => a -> Vector a -> Vector a-{-# INLINE (.>) #-}-s .> v = map (s *) v---- | Vector-scalar division------ >>> fromList [1, 2 :: Double] >/ 2--- fromList [0.5,1.0]-(>/) :: Fractional a => Vector a -> a -> Vector a-{-# INLINE (>/) #-}-v >/ s = map (/ s) v---- | Matrix-scalar division------ >>> fromList [fromList [1, 2], fromList [3, 4 :: Double]] #/ 2--- fromList [fromList [0.5,1.0],fromList [1.5,2.0]]-(#/) :: Fractional a => Vector (Vector a) -> a -> Vector (Vector a)-{-# INLINE (#/) #-}-m #/ s = map (>/ s) m---- | Dot product------ >>> fromList [1, 2] <.> fromList [3, 4 :: Int]--- 11-(<.>) :: Num a => Vector a -> Vector a -> a-{-# INLINE (<.>) #-}-u <.> v = sum $ u >.> v---- | Matrix-vector product------ >>> fromList [fromList [1, 2], fromList [3, 4]] #.> fromList [1, 2 :: Double]--- fromList [5.0,11.0]-(#.>) :: Num a => Vector (Vector a) -> Vector a -> Vector a-{-# INLINE (#.>) #-}-m #.> v = map (<.> v) m---- | Vector-matrix product------ >>> fromList [1, 2 :: Double] <.# fromList [fromList [1, 2], fromList [3, 4]]--- fromList [7.0,10.0]-(<.#) :: Num a => Vector a -> Vector (Vector a) -> Vector a-{-# INLINE (<.#) #-}-v <.# m | any null m = empty-        | otherwise  = (v <.> map head m) `cons` (v <.# map tail m)---- | Matrix transpose------ >>> transpose $ fromList [fromList "ab", fromList "cd"]--- fromList [fromList "ac",fromList "bd"]-transpose :: Vector (Vector a) -> Vector (Vector a)-{-# INLINE transpose #-}-transpose m-  | any null m = empty-  | otherwise  = map head m `cons` transpose (map tail m)
src/Learning/HMM.hs view
@@ -11,7 +11,9 @@  import Prelude hiding (init) import Control.Applicative ((<$>))-import Control.Arrow ((***), first)+import Control.Arrow (first)+import Data.List (elemIndex, genericLength)+import Data.Maybe (fromJust) import Data.Random.Distribution (pdf, rvar) import Data.Random.Distribution.Categorical (Categorical) import qualified Data.Random.Distribution.Categorical as C (@@ -20,15 +22,12 @@ import Data.Random.Distribution.Categorical.Util () import Data.Random.RVar (RVar) import Data.Random.Sample (sample)-import Data.List (genericLength)-import Data.Number.LogFloat (fromLogFloat, logFloat, logFromLogFloat)-import Data.Vector ((!))-import qualified Data.Vector as V (elemIndex, fromList, map, toList, zip)-import qualified Data.Vector.Util.LinearAlgebra as V (transpose)+import qualified Data.Vector as V ((!), elemIndex, fromList, map, toList)+import qualified Data.Vector.Generic as G (convert)+import qualified Data.Vector.Generic.Util.LinearAlgebra as G (transpose)+import qualified Data.Vector.Unboxed as U (fromList, toList) import Learning.HMM.Internal -type LogLikelihood = Double- -- | Parameter set of the hidden Markov model. Direct use of the --   constructor is not recommended. Instead, call 'new' or 'init'. data HMM s o = HMM { states  :: [s] -- ^ Hidden states@@ -89,17 +88,22 @@ -- | @init states outputs@ returns a random variable of the model with --   @states@ and @outputs@, wherein parameters are sampled from uniform --   distributions.-init :: (Ord s, Ord o) => [s] -> [o] -> RVar (HMM s o)-init ss os = do hmm' <- init' (V.fromList ss) (V.fromList os)-                return $ fromHMM' hmm'+init :: (Eq s, Eq o) => [s] -> [o] -> RVar (HMM s o)+init ss os = fromHMM' ss os <$> init' (length ss) (length os)  -- | @model \`withEmission\` xs@ returns a model in which the --   'emissionDist' is updated by using the observed outputs @xs@. The --   'emissionDist' is set to be normalized histograms each of which is --   calculated from a partial set of @xs@ for each state. The partition is --   based on the most likely state path obtained by the Viterbi algorithm.-withEmission :: (Ord s, Ord o) => HMM s o -> [o] -> HMM s o-withEmission model xs = fromHMM' $ withEmission' (toHMM' model) (V.fromList xs)+withEmission :: (Eq s, Eq o) => HMM s o -> [o] -> HMM s o+withEmission model xs = fromHMM' ss os $ withEmission' model' xs'+  where+    ss     = states model+    os     = outputs model+    os'    = V.fromList os+    model' = toHMM' model+    xs'    = U.fromList $ fromJust $ mapM (`V.elemIndex` os') xs  -- | @viterbi model xs@ performs the Viterbi algorithm using the observed --   outputs @xs@, and returns the most likely state path and its log@@ -108,30 +112,36 @@ viterbi model xs =   checkModelIn "viterbi" model `seq`   checkDataIn "viterbi" model xs `seq`-  (V.toList *** logFromLogFloat) $ viterbi' model' xs'+  first (V.toList . V.map (ss V.!) . G.convert) $ viterbi' model' xs'   where+    ss     = V.fromList $ states model+    os'    = V.fromList $ outputs model     model' = toHMM' model-    xs'    = V.fromList xs+    xs'    = U.fromList $ fromJust $ mapM (`V.elemIndex` os') xs --- | @baumWelch model xs@ performs the Baum-Welch algorithm using the---   observed outputs @xs@, and iteratively returns a list of updated---   models and their corresponding log likelihoods.+-- | @baumWelch model xs@ iteratively performs the Baum-Welch algorithm+--   using the observed outputs @xs@, and returns a list of updated models+--   and their corresponding log likelihoods. baumWelch :: (Eq s, Eq o) => HMM s o -> [o] -> [(HMM s o, LogLikelihood)] baumWelch model xs =   checkModelIn "baumWelch" model `seq`   checkDataIn "baumWelch" model xs `seq`-  map (fromHMM' *** logFromLogFloat) $ baumWelch' model' xs'+  map (first $ fromHMM' ss os) $ baumWelch' model' xs'   where+    ss     = states model+    os     = outputs model+    os'    = V.fromList os     model' = toHMM' model-    xs'    = V.fromList xs+    xs'    = U.fromList $ fromJust $ mapM (`V.elemIndex` os') xs  -- | @simulate model t@ generates a Markov process of length @t@ using the --   @model@, and returns its state path and observed outputs. simulate :: HMM s o -> Int -> RVar ([s], [o])-simulate model step | step < 1  = return ([], [])-                    | otherwise = do s0 <- sample $ rvar pi0-                                     x0 <- sample $ rvar $ phi s0-                                     unzip . ((s0, x0) :) <$> sim s0 (step - 1)+simulate model step+  | step < 1  = return ([], [])+  | otherwise = do s0 <- sample $ rvar pi0+                   x0 <- sample $ rvar $ phi s0+                   unzip . ((s0, x0) :) <$> sim s0 (step - 1)   where     sim _ 0 = return []     sim s t = do s' <- sample $ rvar $ w s@@ -164,33 +174,31 @@     err = errorIn fun  -- | Convert 'HMM'' to 'HMM'.-fromHMM' :: (Eq s, Eq o) => HMM' s o -> HMM s o-fromHMM' hmm' = HMM { states           = V.toList ss-                    , outputs          = V.toList os-                    , initialStateDist = C.fromList pi0'-                    , transitionDist   = \s -> case V.elemIndex s ss of-                                                 Nothing -> C.fromList []-                                                 Just i  -> C.fromList $ w' i-                    , emissionDist     = \s -> case V.elemIndex s ss of-                                                 Nothing -> C.fromList []-                                                 Just i  -> C.fromList $ phi' i-                    }+fromHMM' :: (Eq s, Eq o) => [s] -> [o] -> HMM' -> HMM s o+fromHMM' ss os hmm' = HMM { states           = ss+                          , outputs          = os+                          , initialStateDist = C.fromList pi0'+                          , transitionDist   = \s -> case elemIndex s ss of+                                                       Nothing -> C.fromList []+                                                       Just i  -> C.fromList $ w' i+                          , emissionDist     = \s -> case elemIndex s ss of+                                                       Nothing -> C.fromList []+                                                       Just i  -> C.fromList $ phi' i+                          }   where-    ss  = states' hmm'-    os  = outputs' hmm'     pi0 = initialStateDist' hmm'     w   = transitionDist' hmm'-    phi = V.transpose $ emissionDistT' hmm'-    pi0'   = V.toList $ V.map (first fromLogFloat) $ V.zip pi0 ss-    w' i   = V.toList $ V.map (first fromLogFloat) $ V.zip (w ! i) ss-    phi' i = V.toList $ V.map (first fromLogFloat) $ V.zip (phi ! i) os+    phi = G.transpose $ emissionDistT' hmm'+    pi0'   = zip (U.toList pi0) ss+    w' i   = zip (U.toList $ w V.! i) ss+    phi' i = zip (U.toList $ phi V.! i) os  -- | Convert 'HMM' to 'HMM''. The 'initialStateDist'', 'transitionDist'', --   and 'emissionDistT'' are normalized.-toHMM' :: (Eq s, Eq o) => HMM s o -> HMM' s o-toHMM' hmm = HMM' { states'           = V.fromList ss-                  , outputs'          = V.fromList os-                  , initialStateDist' = V.fromList pi0'+toHMM' :: (Eq s, Eq o) => HMM s o -> HMM'+toHMM' hmm = HMM' { nStates'          = length ss+                  , nOutputs'         = length os+                  , initialStateDist' = U.fromList pi0'                   , transitionDist'   = V.fromList w'                   , emissionDistT'    = V.fromList phi'                   }@@ -200,9 +208,9 @@     pi0 = C.normalizeCategoricalPs $ initialStateDist hmm     w   = C.normalizeCategoricalPs . transitionDist hmm     phi = C.normalizeCategoricalPs . emissionDist hmm-    pi0' = [logFloat $ pdf pi0 s | s <- ss]-    w'   = [V.fromList [logFloat $ pdf (w s) s' | s' <- ss] | s <- ss]-    phi' = [V.fromList [logFloat $ pdf (phi s) o | s <- ss] | o <- os]+    pi0' = [pdf pi0 s | s <- ss]+    w'   = [U.fromList [pdf (w s) s' | s' <- ss] | s <- ss]+    phi' = [U.fromList [pdf (phi s) o | s <- ss] | o <- os]  errorIn :: String -> String -> a errorIn fun msg = error $ "Learning.HMM." ++ fun ++ ": " ++ msg
src/Learning/HMM/Internal.hs view
@@ -1,7 +1,6 @@ module Learning.HMM.Internal (     HMM' (..)-  , Likelihood-  , Probability+  , LogLikelihood   , init'   , withEmission'   , viterbi'@@ -12,188 +11,206 @@   ) where  import Control.Applicative ((<$>))+import Control.DeepSeq (NFData, force, rnf) import Control.Monad (forM_, replicateM) import Control.Monad.ST (runST)-import qualified Data.Map.Strict as M (empty, insertWith, findWithDefault)-import Data.Number.LogFloat (LogFloat, logFloat)+import qualified Data.Map.Strict as M (findWithDefault) import Data.Random.RVar (RVar) import Data.Random.Distribution.Simplex (stdSimplex)-import Data.Random.Distribution.Uniform.Util ()-import Data.Vector (Vector, (!)) import qualified Data.Vector as V (-    filter, foldl', foldl1', freeze, fromList, last, length, map, maximum-  , maxIndex, replicate, sum, tail, zip, zipWith, zipWith3, zipWith4+    Vector, (!), filter, foldl1', freeze, fromList, generate, map , tail+  , zip, zipWith, zipWith3   )+import qualified Data.Vector.Generic as G (convert)+import qualified Data.Vector.Generic.Util as G (frequencies)+import Data.Vector.Generic.Util.LinearAlgebra (+    (>+>), (>.>), (>/>), (#+#), (.>), (>/), (#.>), (<.#)+  )+import qualified Data.Vector.Generic.Util.LinearAlgebra as G (transpose) import qualified Data.Vector.Mutable as MV (new, read, write)-import qualified Data.Vector.Util as V (unsafeElemIndex)-import Data.Vector.Util.LinearAlgebra (-    (>+>), (>.>), (>/>), (#+#), (.>), (>/), (#/), (<.>), (#.>), (<.#)+import qualified Data.Vector.Unboxed as U (+    Vector, (!), freeze, fromList, generate, length, map, maxIndex, maximum+  , replicate, sum, tail, zip   )-import qualified Data.Vector.Util.LinearAlgebra as V (transpose)+import qualified Data.Vector.Unboxed.Mutable as MU (new, read, write) -type Likelihood  = LogFloat-type Probability = LogFloat+type LogLikelihood = Double --- | More efficient data structure of the 'HMM' model. This should be---   only used internally. The 'emissionDistT'' is a transposed matrix in---   order to simplify the calculation.-data HMM' s o = HMM' { states'           :: Vector s-                     , outputs'          :: Vector o-                     , initialStateDist' :: Vector Probability-                     , transitionDist'   :: Vector (Vector Probability)-                     , emissionDistT'    :: Vector (Vector Probability)-                     }+-- | More efficient data structure of the 'HMM' model. The 'states' and+--   'outputs' in 'HMM' are represented by their indices. The+--   'initialStateDist', 'transitionDist', and 'emissionDist' are+--   represented by matrices. The 'emissionDistT'' is a transposed matrix+--   in order to simplify the calculation.+data HMM' = HMM' { nStates'          :: Int -- ^ Number of states+                 , nOutputs'         :: Int -- ^ Number of outputs+                 , initialStateDist' :: U.Vector Double+                 , transitionDist'   :: V.Vector (U.Vector Double)+                 , emissionDistT'    :: V.Vector (U.Vector Double)+                 } -init' :: Vector s -> Vector o -> RVar (HMM' s o)-init' ss os = do-  let n = V.length ss-      m = V.length os-  pi0 <- V.fromList <$> stdSimplex (n-1)-  w   <- V.fromList <$> replicateM n (V.fromList <$> stdSimplex (n-1))-  phi <- V.fromList <$> replicateM n (V.fromList <$> stdSimplex (m-1))-  return HMM' { states'           = ss-              , outputs'          = os+instance NFData HMM' where+  rnf hmm' = rnf n `seq` rnf m `seq` rnf pi0 `seq` rnf w `seq` rnf phi'+    where+      n    = nStates' hmm'+      m    = nOutputs' hmm'+      pi0  = initialStateDist' hmm'+      w    = transitionDist' hmm'+      phi' = emissionDistT' hmm'++init' :: Int -> Int -> RVar HMM'+init' n m = do+  pi0 <- U.fromList <$> stdSimplex (n-1)+  w   <- V.fromList <$> replicateM n (U.fromList <$> stdSimplex (n-1))+  phi <- V.fromList <$> replicateM n (U.fromList <$> stdSimplex (m-1))+  return HMM' { nStates'          = n+              , nOutputs'         = m               , initialStateDist' = pi0               , transitionDist'   = w-              , emissionDistT'    = V.transpose phi+              , emissionDistT'    = G.transpose phi               } -withEmission' :: (Ord s, Ord o) => HMM' s o -> Vector o -> HMM' s o+withEmission' :: HMM' -> U.Vector Int -> HMM' withEmission' model xs = model { emissionDistT' = phi' }   where-    ss = states' model-    os = outputs' model-    (path, _) = viterbi' model xs-    mp    = V.foldl' (\m k -> M.insertWith (+) k 1 m) M.empty $ V.zip path xs-    hists = V.map (\s -> V.map (\o -> M.findWithDefault 0 (s, o) mp) os) ss-    phi'  = V.transpose $ V.map (\h -> h >/ V.sum h) hists+    ss   = V.generate (nStates' model) id+    os   = U.generate (nOutputs' model) id+    phi' = let (path, _) = viterbi' model xs+               freqs     = G.frequencies $ U.zip path xs+               hists     = V.map (\s -> U.map (\o ->+                                 M.findWithDefault 0 (s, o) freqs) os) ss+           in V.map (\f -> f >/ U.sum f) hists -viterbi' :: Eq o => HMM' s o -> Vector o -> (Vector s, Likelihood)-viterbi' model xs = (path, likelihood)+viterbi' :: HMM' -> U.Vector Int -> (U.Vector Int, LogLikelihood)+viterbi' model xs = (path, logL)   where-    -- The following procedure is based on-    -- http://ibisforest.org/index.php?cmd=read&page=Viterbi%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0&word=Viterbi-    path = V.map (ss !) $ runST $ do-      ix <- MV.new n-      ix `MV.write` (n-1) $ V.maxIndex $ deltas ! (n-1)-      forM_ (reverse [0..(n-2)]) $ \i -> do-        j <- ix `MV.read` (i+1)-        ix `MV.write` i $ psis ! (i+1) ! j-      V.freeze ix-      where-        ss = states' model-    likelihood = V.maximum $ deltas ! (n-1)+    n = U.length xs -    deltas :: Vector (Vector Probability)-    psis   :: Vector (Vector Int)+    -- First, we calculate the value function and the state maximizing it+    -- for each time.+    deltas :: V.Vector (U.Vector Double)+    psis   :: V.Vector (U.Vector Int)     (deltas, psis) = runST $ do       ds <- MV.new n       ps <- MV.new n-      ds `MV.write` 0 $ (phi' ! x 0) >.> pi0-      ps `MV.write` 0 $ V.replicate k (0 :: Int)-      forM_ [1..(n-1)] $ \i -> do-        d <- ds `MV.read` (i-1)-        let dws = V.map (d >.>) w'-        ds `MV.write` i $ phi' ! x i >.> V.map V.maximum dws-        ps `MV.write` i $ V.map V.maxIndex dws+      MV.write ds 0 $ U.map log (phi' V.! (xs U.! 0)) >+> U.map log pi0+      MV.write ps 0 $ U.replicate k 0+      forM_ [1..(n-1)] $ \t -> do+        d <- MV.read ds (t-1)+        let dws = V.map (\wj -> d >+> U.map log wj) w'+        MV.write ds t $ U.map log (phi' V.! (xs U.! t)) >+> G.convert (V.map U.maximum dws)+        MV.write ps t $ G.convert (V.map U.maxIndex dws)       ds' <- V.freeze ds       ps' <- V.freeze ps       return (ds', ps')       where-        k   = V.length $ states' model-        x i = let os  = outputs' model-                  xs' = V.map (`V.unsafeElemIndex` os) xs-              in xs' ! i+        k    = nStates' model         pi0  = initialStateDist' model-        w'   = V.transpose $ transitionDist' model+        w'   = G.transpose $ transitionDist' model         phi' = emissionDistT' model -    -- Here we assumed that-    n = V.length xs+    -- The most likely path and corresponding log likelihood is as follows.+    path = runST $ do+      ix <- MU.new n+      MU.write ix (n-1) $ U.maxIndex (deltas V.! (n-1))+      forM_ (reverse [0..(n-2)]) $ \t -> do+        i <- MU.read ix (t+1)+        MU.write ix t $ psis V.! (t+1) U.! i+      U.freeze ix+    logL = U.maximum $ deltas V.! (n-1) -baumWelch' :: (Eq s, Eq o) => HMM' s o -> Vector o -> [(HMM' s o, Likelihood)]-baumWelch' model xs = zip ms $ tail ells+baumWelch' :: HMM' -> U.Vector Int -> [(HMM', LogLikelihood)]+baumWelch' model xs = zip models (tail logLs)   where-    (ms, ells) = unzip $ iterate ((`baumWelch1'` xs) . fst) (model, undefined)+    n = U.length xs+    step (m, _)     = baumWelch1' m n xs+    (models, logLs) = unzip $ iterate step (model, undefined)  -- | Perform one step of the Baum-Welch algorithm and return the updated --   model and the likelihood of the old model.-baumWelch1' :: (Eq s, Eq o) => HMM' s o -> Vector o -> (HMM' s o, Likelihood)-baumWelch1' model xs = (model', likelihood)+baumWelch1' :: HMM' -> Int -> U.Vector Int -> (HMM', LogLikelihood)+baumWelch1' model n xs = force (model', logL)   where-    model' = model { initialStateDist' = pi0-                   , transitionDist'   = w-                   , emissionDistT'    = phi'-                   }-    likelihood = V.last ells--    -- First, we calculate the alpha and beta values using the+    -- First, we calculate the alpha, beta, and scaling values using the     -- forward-backward algorithm.-    alphas = forward' model xs-    betas  = backward' model xs+    (alphas, cs) = forward' model n xs+    betas        = backward' model n xs cs -    -- Then, we obtain the likelihoods for each time. This should be-    -- constant over time.-    ells = V.zipWith (<.>) alphas betas+    -- Based on the alpha, beta, and scaling values, we calculate the+    -- posterior distribution, i.e., gamma and xi values.+    (gammas, xis) = posterior' model n xs alphas betas cs -    -- Based on the alpha, beta, and likelihood values, we calculate the-    -- gamma and xi values.-    gammas = V.zipWith3 (\a b l -> a >.> b >/ l) alphas betas ells-    xis    = V.zipWith4 (\a b l x -> let w1 = V.zipWith (.>) a w0-                                         w2 = V.map (phi0 ! x >.> b >.>) w1-                                     in w2 #/ l)-               alphas (V.tail betas) (V.tail ells) (V.tail xs')+    -- Using the gamma and xi values, we obtain the optimal initial state+    -- probability vector, transition probability matrix, and emission+    -- probability matrix.+    pi0  = gammas V.! 0+    w    = let ds = V.foldl1' (#+#) xis -- denominators+               ns = V.map U.sum ds      -- numerators+           in V.zipWith (>/) ds ns+    phi' = let gs' o = V.map snd $ V.filter ((== o) . fst) $ V.zip (G.convert xs) gammas+               ds    = V.foldl1' (>+>) . gs'  -- denominators+               ns    = V.foldl1' (>+>) gammas -- numerators+           in V.map (\o -> ds o >/> ns) os       where-        xs'  = V.map (`V.unsafeElemIndex` os) xs-        w0   = transitionDist' model-        phi0 = emissionDistT' model--    -- Using the gamma and xi values, we finally obtain the optimal initial-    -- state probability vector, transition probability matrix, and-    -- emission probability matrix.-    pi0  = let gs = gammas ! 0-           in gs >/ V.sum gs-    w    = let ws = V.foldl1' (#+#) xis-               zs = V.map V.sum ws-           in V.zipWith (>/) ws zs-    phi' = let gs' o = V.map snd $ V.filter ((== o) . fst) $ V.zip xs gammas-               phis  = V.foldl1' (>+>) . gs'-               zs    = V.foldl1' (>+>) gammas-           in V.map (\o -> phis o >/> zs) os+        os = V.generate (nOutputs' model) id -    -- Here we assumed that-    os = outputs' model+    -- We finally obtain the new model and the likelihood for the old model.+    model' = model { initialStateDist' = pi0+                   , transitionDist'   = w+                   , emissionDistT'    = phi'+                   }+    logL = - (U.sum $ U.map log cs) -forward' :: Eq o => HMM' s o -> Vector o -> Vector (Vector Probability)-forward' model xs = runST $ do-  v <- MV.new n-  v `MV.write` 0 $ (phi' ! x 0) >.> pi0-  forM_ [1..(n-1)] $ \i -> do-    a <- v `MV.read` (i-1)-    v `MV.write` i $ (phi' ! x i) >.> (a <.# w)-  V.freeze v+-- | Return alphas and scaling variables.+forward' :: HMM' -> Int -> U.Vector Int -> (V.Vector (U.Vector Double), U.Vector Double)+{-# INLINE forward' #-}+forward' model n xs = runST $ do+  as <- MV.new n+  cs <- MU.new n+  let a0 = (phi' V.! (xs U.! 0)) >.> pi0+      c0 = 1 / U.sum a0+  MV.write as 0 (c0 .> a0)+  MU.write cs 0 c0+  forM_ [1..(n-1)] $ \t -> do+    a <- MV.read as (t-1)+    let a' = (phi' V.! (xs U.! t)) >.> (a <.# w)+        c' = 1 / U.sum a'+    MV.write as t (c' .> a')+    MU.write cs t c'+  as' <- V.freeze as+  cs' <- U.freeze cs+  return (as', cs')   where-    n   = V.length xs-    x i = let os  = outputs' model-              xs' = V.map (`V.unsafeElemIndex` os) xs-          in xs' ! i     pi0  = initialStateDist' model     w    = transitionDist' model     phi' = emissionDistT' model -backward' :: Eq o => HMM' s o -> Vector o -> Vector (Vector Probability)-backward' model xs = runST $ do-  v <- MV.new n-  v `MV.write` (n-1) $ V.replicate k $ logFloat (1 :: Double)-  forM_ (reverse [0..(n-2)]) $ \i -> do-    b <- v `MV.read` (i+1)-    v `MV.write` i $ w #.> ((phi' ! x (i+1)) >.> b)-  V.freeze v+-- | Return betas using scaling variables.+backward' :: HMM' -> Int -> U.Vector Int -> U.Vector Double -> V.Vector (U.Vector Double)+{-# INLINE backward' #-}+backward' model n xs cs = runST $ do+  bs <- MV.new n+  let bE = U.replicate k 1+      cE = cs U.! (n-1)+  MV.write bs (n-1) $ cE .> bE+  forM_ (reverse [0..(n-2)]) $ \t -> do+    b <- MV.read bs (t+1)+    let b' = w #.> ((phi' V.! (xs U.! (t+1))) >.> b)+        c' = cs U.! t+    MV.write bs t $ c' .> b'+  V.freeze bs   where-    n   = V.length xs-    k   = V.length $ states' model-    x i = let os  = outputs' model-              xs' = V.map (`V.unsafeElemIndex` os) xs-          in xs' ! i+    k    = nStates' model+    w    = transitionDist' model+    phi' = emissionDistT' model++-- | Return the posterior distribution.+posterior' :: HMM' -> Int -> U.Vector Int -> V.Vector (U.Vector Double) -> V.Vector (U.Vector Double) -> U.Vector Double -> (V.Vector (U.Vector Double), V.Vector (V.Vector (U.Vector Double)))+{-# INLINE posterior' #-}+posterior' model _ xs alphas betas cs = (gammas, xis)+  where+    gammas = V.zipWith3 (\a b c -> a >.> b >/ c) alphas betas (G.convert cs)+    xis    = V.zipWith3 (\a b x -> let w' = V.zipWith (.>) (G.convert a) w+                                   in V.map ((phi' V.! x) >.> b >.>) w')+               alphas (V.tail betas) (G.convert $ U.tail xs)     w    = transitionDist' model     phi' = emissionDistT' model
tests/doctests.hs view
@@ -4,6 +4,6 @@  main :: IO () main = doctest [ "-isrc"-               , "src/Data/Vector/Util/LinearAlgebra.hs"+               , "src/Data/Vector/Generic/Util/LinearAlgebra.hs"                , "src/Learning/HMM.hs"                ]