packages feed

hmm-hmatrix (empty) → 0.0

raw patch · 15 files changed

+1941/−0 lines, 15 filesdep +arraydep +basedep +containerssetup-changed

Dependencies added: array, base, containers, explicit-exception, hmatrix, lazy-csv, non-empty, random, semigroups, transformers, utility-ht

Files

+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2015, Henning Thielemann++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 Henning Thielemann 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.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ hmm-hmatrix.cabal view
@@ -0,0 +1,79 @@+Name:                hmm-hmatrix+Version:             0.0+Synopsis:            Hidden Markov Models using HMatrix primitives+Description:+  Hidden Markov Models implemented using HMatrix data types and operations.+  <http://en.wikipedia.org/wiki/Hidden_Markov_Model>+  .+  It implements:+  .+  * generation of samples of emission sequences,+  .+  * computation of the likelihood of an observed sequence of emissions,+  .+  * construction of most likely state sequence+    that produces an observed sequence of emissions,+  .+  * supervised and unsupervised training of the model by Baum-Welch algorithm.+  .+  It supports any kind of emission distribution,+  where discrete and multivariate Gaussian distributions+  are implemented as examples.+  .+  For an introduction please refer to the examples:+  .+  * "Math.HiddenMarkovModel.Example.TrafficLight"+  .+  * "Math.HiddenMarkovModel.Example.SineWave"+  .+  * "Math.HiddenMarkovModel.Example.Circle"+  .+  An alternative package without foreign calls is @hmm@.+Homepage:            http://code.haskell.org/~thielema/hmm-hmatrix+License:             BSD3+License-File:        LICENSE+Author:              Henning Thielemann+Maintainer:          haskell@henning-thielemann.de+Category:            Math+Build-Type:          Simple+Cabal-Version:       >=1.10++Source-Repository this+  Tag:         0.0+  Type:        darcs+  Location:    http://code.haskell.org/~thielema/hmm-hmatrix++Source-Repository head+  Type:        darcs+  Location:    http://code.haskell.org/~thielema/hmm-hmatrix++Library+  Exposed-Modules:+    Math.HiddenMarkovModel+    Math.HiddenMarkovModel.Named+    Math.HiddenMarkovModel.Distribution+    Math.HiddenMarkovModel.Pattern+    Math.HiddenMarkovModel.Example.TrafficLight+    Math.HiddenMarkovModel.Example.SineWave+    Math.HiddenMarkovModel.Example.Circle+  Other-Modules:+    Math.HiddenMarkovModel.Normalized+    Math.HiddenMarkovModel.Private+    Math.HiddenMarkovModel.Utility+    Math.HiddenMarkovModel.CSV+    Math.HiddenMarkovModel.Test+  Build-Depends:+    hmatrix >=0.15 && <0.16,+    explicit-exception >=0.1.7 && <0.2,+    lazy-csv >=0.5 && <0.6,+    random >=1.0 && <1.1,+    transformers >= 0.2 && <0.5,+    non-empty >=0.2.1 && <0.3,+    semigroups >=0.8.4.1 && <0.17,+    containers >=0.4.2 && <0.6,+    array >=0.4 && <0.6,+    utility-ht >=0.0.10 && <0.1,+    base >=4.5 && <4.8+  Hs-Source-Dirs:      src+  Default-Language:    Haskell2010+  GHC-Options:         -Wall
+ src/Math/HiddenMarkovModel.hs view
@@ -0,0 +1,178 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+module Math.HiddenMarkovModel (+   T(..), Distr.State, state,+   Discrete, DiscreteTrained,+   Gaussian, GaussianTrained,+   uniform,+   generate,+   Normalized.logLikelihood,+   Normalized.reveal,++   Trained(..),+   trainSupervised,+   Normalized.trainUnsupervised,+   mergeTrained, finishTraining, trainMany,+   deviation,++   toCSV,+   fromCSV,+   ) where++import qualified Math.HiddenMarkovModel.Distribution as Distr+import qualified Math.HiddenMarkovModel.Normalized as Normalized+import qualified Math.HiddenMarkovModel.CSV as HMMCSV+import Math.HiddenMarkovModel.Private+          (T(..), Trained(..), mergeTrained, toCells, parseCSV)+import Math.HiddenMarkovModel.Distribution (State(State))+import Math.HiddenMarkovModel.Utility+          (randomItemProp, normalizeProb, attachOnes)++import qualified Numeric.LinearAlgebra.Algorithms as Algo+import qualified Numeric.Container as NC+import qualified Data.Packed.Matrix as Matrix+import qualified Data.Packed.Vector as Vector+import Data.Packed.Matrix (Matrix)+import Data.Packed.Vector (Vector)++import qualified Text.CSV.Lazy.String as CSV++import qualified System.Random as Rnd++import qualified Control.Monad.Exception.Synchronous as ME+import qualified Control.Monad.Trans.State as MS+import qualified Control.Monad.HT as Monad++import qualified Data.NonEmpty as NonEmpty+import qualified Data.Array as Array+import Data.Foldable (Foldable)+import Data.Array (accumArray)++++state :: Int -> State+state = State+++type DiscreteTrained prob symbol = Trained (Distr.DiscreteTrained prob symbol) prob+type Discrete prob symbol = T (Distr.Discrete prob symbol) prob++type GaussianTrained a = Trained (Distr.GaussianTrained a) a+type Gaussian a = T (Distr.Gaussian a) a+++{- |+Create a model with uniform probabilities+for initial vector and transition matrix+given a distribution for the emissions.+You can use this as a starting point for 'Normalized.trainUnsupervised'.+-}+uniform ::+   (Distr.Info distr, Distr.Probability distr ~ prob) =>+   distr -> T distr prob+uniform distr =+   let n = Distr.numberOfStates distr+       c = recip $ fromIntegral n+   in  Cons {+          initial = NC.constant c n,+          transition = NC.konst c (n,n),+          distribution = distr+       }+++generate ::+   (Rnd.RandomGen g, Ord prob, Rnd.Random prob,+    Distr.Generate distr, Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   T distr prob -> g -> [emission]+generate hmm =+   MS.evalState $+   flip MS.evalStateT (initial hmm) $+   Monad.repeat $ MS.StateT $ \v0 -> do+      s <- randomItemProp $ zip [0..] (Vector.toList v0)+      x <- Distr.generate (distribution hmm) (State s)+      return (x, takeColumn s $ transition hmm)++takeColumn :: (Matrix.Element a) => Int -> Matrix a -> Vector a+takeColumn n  =  Matrix.flatten . Matrix.extractRows [n] . Matrix.trans++++{- |+Contribute a manually labeled emission sequence to a HMM training.+-}+trainSupervised ::+   (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   Int -> NonEmpty.T [] (State, emission) -> Trained tdistr prob+trainSupervised n xs =+   let getState (State s, _x) = s+   in  Trained {+          trainedInitial = NC.assoc n 0 [(getState (NonEmpty.head xs), 1)],+          trainedTransition =+             Matrix.trans $ NC.accum (NC.konst 0 (n,n)) (+) $+             attachOnes $ NonEmpty.mapAdjacent (,) $ fmap getState xs,+          trainedDistribution =+             Distr.accumulateEmissions $ map attachOnes $ Array.elems $+             accumArray (flip (:)) [] (State 0, State (n-1)) $ NonEmpty.flatten xs+       }++finishTraining ::+   (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   Trained tdistr prob -> T distr prob+finishTraining hmm =+   Cons {+      initial = normalizeProb $ trainedInitial hmm,+      transition =+         Matrix.fromColumns $ map normalizeProb $+         Matrix.toColumns $ trainedTransition hmm,+      distribution = Distr.normalize $ trainedDistribution hmm+   }++trainMany ::+   (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,+    Distr.Probability distr ~ prob,+    Foldable f) =>+   (trainingData -> Trained tdistr prob) ->+   NonEmpty.T f trainingData -> T distr prob+trainMany train =+   finishTraining . NonEmpty.foldl1Map mergeTrained train++++++{- |+Compute maximum deviation between initial and transition probabilities.+You can use this as abort criterion for unsupervised training.+We omit computation of differences between the emission probabilities.+This simplifies matters a lot and+should suffice for defining an abort criterion.+-}+deviation ::+   (Algo.Field prob, Ord prob) => T distr prob -> T distr prob -> prob+deviation hmm0 hmm1 =+   deviationVec (initial hmm0) (initial hmm1)+   `max`+   deviationVec (transition hmm0) (transition hmm1)++deviationVec ::+   (Ord a, NC.Container c a) =>+   c a -> c a -> a+deviationVec x y =+   let d = NC.sub x y+   in  NC.maxElement d `max` negate (NC.minElement d)+++toCSV ::+   (Distr.CSV distr, Algo.Field prob, Show prob) =>+   T distr prob -> String+toCSV hmm =+   CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $+   toCells hmm++fromCSV ::+   (Distr.CSV distr, Algo.Field prob, Read prob) =>+   String -> ME.Exceptional String (T distr prob)+fromCSV =+   MS.evalStateT parseCSV . map HMMCSV.fixShortRow . CSV.parseCSV
+ src/Math/HiddenMarkovModel/CSV.hs view
@@ -0,0 +1,151 @@+module Math.HiddenMarkovModel.CSV where++import qualified Numeric.LinearAlgebra.Algorithms as Algo+import qualified Data.Packed.Matrix as Matrix+import qualified Data.Packed.Vector as Vector+import Data.Packed.Matrix (Matrix)+import Data.Packed.Vector (Vector)++import qualified Text.CSV.Lazy.String as CSV+import Text.Read.HT (maybeRead)+import Text.Printf (printf)++import qualified Control.Monad.Exception.Synchronous as ME+import qualified Control.Monad.Trans.Class as MT+import qualified Control.Monad.Trans.State as MS+import Control.Monad.Exception.Synchronous (Exceptional)+import Control.Monad (liftM2, replicateM, unless)++import qualified Data.List.HT as ListHT+++cellsFromVector ::+   (Show a, Algo.Field a) =>+   Vector a -> [String]+cellsFromVector = map show . Vector.toList++cellsFromMatrix ::+   (Show a, Matrix.Element a) =>+   Matrix.Matrix a -> [[String]]+cellsFromMatrix = map (map show) . Matrix.toLists++padTable :: a -> [[a]] -> [[a]]+padTable x xs =+   let width = maximum (map length xs)+   in  map (ListHT.padRight x width) xs+++type CSVParser = MS.StateT CSV.CSVResult (Exceptional String)++assert :: Bool -> String -> CSVParser ()+assert cond msg =+   unless cond $ MT.lift $ ME.throw msg++retrieveShortRow :: CSV.CSVError -> Maybe CSV.CSVRow+retrieveShortRow err =+   case err of+      CSV.IncorrectRow {CSV.csvFields = row} -> Just row+      _ -> Nothing++fixShortRow ::+   Either [CSV.CSVError] CSV.CSVRow -> Either [CSV.CSVError] CSV.CSVRow+fixShortRow erow =+   case erow of+      Left errs ->+         case ListHT.partitionMaybe retrieveShortRow errs of+            ([row], []) -> Right row+            _ -> Left errs+      _ -> erow++maybeGetRow :: CSVParser (Maybe CSV.CSVRow)+maybeGetRow = do+   csv0 <- MS.get+   case csv0 of+      [] -> return Nothing+      item : csv1 -> do+         MS.put csv1+         case item of+            Right row -> return (Just row)+            Left errors ->+               MT.lift $ ME.throw $ unlines $ map CSV.ppCSVError errors++getRow :: CSVParser CSV.CSVRow+getRow =+   MT.lift . ME.fromMaybe "unexpected end of file" =<< maybeGetRow++checkEmptyRow :: CSV.CSVRow -> Exceptional String ()+checkEmptyRow row =+   case filter (not . null . CSV.csvFieldContent) row of+      [] -> return ()+      cell:_ -> ME.throw $ printf "%d: expected empty row" (CSV.csvRowNum cell)++skipEmptyRow :: CSVParser ()+skipEmptyRow  =  MT.lift . checkEmptyRow =<< getRow++manySepUntilEnd :: CSVParser a -> CSVParser [a]+manySepUntilEnd p =+   let go = liftM2 (:) p $ do+          mrow <- maybeGetRow+          case mrow of+             Nothing -> return []+             Just row -> do+                MT.lift $ checkEmptyRow row+                go+   in  go++manyRowsUntilEnd :: (CSV.CSVRow -> CSVParser a) -> CSVParser [a]+manyRowsUntilEnd p =+   let go = do+          mrow <- maybeGetRow+          case mrow of+             Nothing -> return []+             Just row -> liftM2 (:) (p row) go+   in  go++parseVectorCells ::+   (Read a, Algo.Field a) =>+   CSVParser (Vector a)+parseVectorCells =+   parseVectorFields =<< getRow++parseVectorFields ::+   (Read a, Algo.Field a) =>+   CSV.CSVRow -> CSVParser (Vector a)+parseVectorFields =+   MT.lift . fmap Vector.fromList . mapM parseNumberCell .+   ListHT.dropWhileRev (null . CSV.csvFieldContent)++parseNonEmptyVectorCells ::+   (Read a, Algo.Field a) =>+   CSVParser (Vector a)+parseNonEmptyVectorCells = do+   v <- parseVectorCells+   assert (Vector.dim v > 0) "no data for vector"+   return v++cellContent :: CSV.CSVField -> Exceptional String String+cellContent field =+   case field of+      CSV.CSVFieldError {} -> ME.throw $ CSV.ppCSVField field+      CSV.CSVField { CSV.csvFieldContent = str } -> return str++parseNumberCell :: (Read a) => CSV.CSVField -> Exceptional String a+parseNumberCell field = do+   str <- cellContent field+   ME.fromMaybe (printf "field content \"%s\" is not a number" str) $+      maybeRead str++parseSquareMatrixCells ::+   (Read a, Algo.Field a) =>+   Int -> CSVParser (Matrix a)+parseSquareMatrixCells n = do+   rows <- replicateM n parseVectorCells+   assert (not $ null rows) "no rows"+   assert (all ((n==) . Vector.dim) rows) "inconsistent matrix dimensions"+   return $ Matrix.fromRows rows++parseStringList ::+   CSV.CSVRow -> CSVParser [String]+parseStringList =+   MT.lift . mapM cellContent .+   ListHT.dropWhileRev (null . CSV.csvFieldContent)
+ src/Math/HiddenMarkovModel/Distribution.hs view
@@ -0,0 +1,327 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+module Math.HiddenMarkovModel.Distribution (+   State(..),+   Emission, Probability, Trained,+   Info(..), Generate(..), EmissionProb(..), Estimate(..),++   Discrete(..), DiscreteTrained(..),+   Gaussian(..), GaussianTrained(..), gaussian,++   CSV(..), HMMCSV.CSVParser, CSVSymbol(..),+   ) where++import qualified Math.HiddenMarkovModel.CSV as HMMCSV+import Math.HiddenMarkovModel.Utility (randomItemProp, normalizeProb)++import qualified Numeric.LinearAlgebra.Algorithms as Algo+import qualified Numeric.Container as NC+import qualified Data.Packed.Matrix as Matrix+import qualified Data.Packed.Vector as Vector+import Numeric.Container ((<>))+import Data.Packed.Matrix (Matrix)+import Data.Packed.Vector (Vector)++import qualified System.Random as Rnd++import qualified Text.CSV.Lazy.String as CSV+import Text.Read.HT (maybeRead)+import Text.Printf (printf)++import qualified Control.Monad.Exception.Synchronous as ME+import qualified Control.Monad.Trans.Class as MT+import qualified Control.Monad.Trans.State as MS+import Control.Monad (liftM2)++import qualified Data.NonEmpty as NonEmpty+import qualified Data.Foldable as Fold+import qualified Data.Map as Map+import qualified Data.Set as Set+import qualified Data.Array as Array+import qualified Data.List as List+import Data.Foldable (foldMap)+import Data.Tuple.HT (mapFst)+import Data.Array (Array, Ix, listArray, (!))+import Data.Map (Map)+import Data.Maybe (listToMaybe)+++newtype State = State Int+   deriving (Eq, Ord, Show, Read, Ix)++instance Enum State where+   toEnum = State+   fromEnum (State n) = n+++type family Probability distr+type family Emission distr+type family Trained distr+++class+   (NC.Container Vector (Probability distr), NC.Product (Probability distr)) =>+      Info distr where+   numberOfStates :: distr -> Int++class+   (NC.Container Vector (Probability distr), NC.Product (Probability distr)) =>+      Generate distr where+   generate ::+      (Rnd.RandomGen g, Probability distr ~ prob, Emission distr ~ emission) =>+      distr -> State -> MS.State g emission++class+   (NC.Container Vector (Probability distr), NC.Product (Probability distr)) =>+      EmissionProb distr where+   emissionProb :: distr -> Emission distr -> Vector (Probability distr)++class+   (EmissionProb (Distribution tdistr),+    Trained (Distribution tdistr) ~ tdistr) =>+      Estimate tdistr where+   type Distribution tdistr+   accumulateEmissions ::+      (Distribution tdistr ~ distr, Probability distr ~ prob) =>+      [[(Emission distr, prob)]] -> tdistr+   -- could as well be in Semigroup class+   combine :: tdistr -> tdistr -> tdistr+   normalize :: (Distribution tdistr ~ distr) => tdistr -> distr++++newtype Discrete prob symbol = Discrete (Map symbol (Vector prob))+   deriving (Show)++newtype DiscreteTrained prob symbol = DiscreteTrained (Map symbol (Vector prob))+   deriving (Show)++type instance Probability (Discrete prob symbol) = prob+type instance Emission (Discrete prob symbol) = symbol++type instance Trained (Discrete prob symbol) = DiscreteTrained prob symbol++instance+   (NC.Container Vector prob, NC.Product prob, Ord symbol) =>+      Info (Discrete prob symbol) where+   numberOfStates (Discrete m) = Vector.dim $ snd $ Map.findMin m++instance+   (NC.Container Vector prob, NC.Product prob, Ord symbol,+    Ord prob, Rnd.Random prob) =>+      Generate (Discrete prob symbol) where+   generate (Discrete m) (State state) =+      randomItemProp $ Map.toAscList $ fmap (flip NC.atIndex state) m++instance+   (NC.Container Vector prob, NC.Product prob, Ord symbol) =>+      EmissionProb (Discrete prob symbol) where+   emissionProb (Discrete m) =+      mapLookup "emitDiscrete: unknown emission symbol" m++instance+   (NC.Container Vector prob, NC.Product prob, Ord symbol) =>+      Estimate (DiscreteTrained prob symbol) where+   type Distribution (DiscreteTrained prob symbol) = Discrete prob symbol+   accumulateEmissions grouped =+      let set = Set.toAscList $ foldMap (Set.fromList . map fst) grouped+          emi = Map.fromAscList $ zip set [0..]+      in  DiscreteTrained $ Map.fromAscList $ zip set $+          transposeVectorList $+          map+             (NC.accum (NC.konst 0 (length set)) (+) .+              map (mapFst+                 (mapLookup "estimateDiscrete: unknown emission symbol" emi)))+             grouped+   combine (DiscreteTrained distr0) (DiscreteTrained distr1) =+      DiscreteTrained $ Map.unionWith NC.add distr0 distr1+   normalize (DiscreteTrained distr) =+      Discrete $ Map.fromAscList $ zip (Map.keys distr) $+      transposeVectorList $ map normalizeProb $+      transposeVectorList $ Map.elems distr++transposeVectorList :: (Matrix.Element a) => [Vector a] -> [Vector a]+transposeVectorList = Matrix.toRows . Matrix.fromColumns++mapLookup :: (Ord k) => String -> Map.Map k a -> k -> a+mapLookup msg dict x =+   Map.findWithDefault (error msg) x dict+++newtype Gaussian a = Gaussian (Array State (Vector a, Matrix a, a))+   deriving (Show)++newtype GaussianTrained a = GaussianTrained (Map State (Vector a, Matrix a, a))+   deriving (Show)++type instance Probability (Gaussian a) = a+type instance Emission (Gaussian a) = Vector a++type instance Trained (Gaussian a) = GaussianTrained a++instance (Algo.Field a) => Info (Gaussian a) where+   numberOfStates (Gaussian params) = Array.rangeSize $ Array.bounds params++instance (Algo.Field a, Ord a, Rnd.Random a) => Generate (Gaussian a) where+   generate (Gaussian allParams) state = do+      let (center, covarianceChol, _c) = allParams ! state+      seed <- MS.state Rnd.random+      return $+         NC.add center $+         NC.cmap realToFrac+               (NC.randomVector seed NC.Gaussian (Vector.dim center))+            <> covarianceChol++instance (Algo.Field a) => EmissionProb (Gaussian a) where+   emissionProb (Gaussian allParams) =+      let cholSolve m x =+             Matrix.flatten $ Algo.cholSolve m $ Matrix.asColumn x+          prob x (center, covarianceChol, c) =+             let x0 = NC.sub x center+             in  c * exp ((-1/2) * NC.dot x0 (cholSolve covarianceChol x0))+      in  \x -> Vector.fromList $ map (prob x) $ Array.elems allParams++instance (Algo.Field a) => Estimate (GaussianTrained a) where+   type Distribution (GaussianTrained a) = Gaussian a+   accumulateEmissions =+      let params xs =+             let center =+                    NonEmpty.foldl1Map NC.add (uncurry $ flip NC.scale) xs+                 covariance =+                    NonEmpty.foldl1Map NC.add (\(x,c) -> NC.scale c $ NC.outer x x) xs+             in  (center, covariance, Fold.sum $ fmap snd xs)+      in  GaussianTrained . fmap params . Map.mapMaybe NonEmpty.fetch .+          Map.fromList . zip [State 0 ..]+   combine (GaussianTrained distr0) (GaussianTrained distr1) =+      let comb (center0, covariance0, weight0)+               (center1, covariance1, weight1) =+             (NC.add center0 center1,+              NC.add covariance0 covariance1,+              weight0 + weight1)+      in  GaussianTrained $ Map.unionWith comb distr0 distr1+   {-+     Sum_i (xi-mi) * (xi-mi)^T+   = Sum_i xi*xi^T + Sum_i mi*mi^T - Sum_i xi*mi^T - Sum_i mi*xi^T+   = Sum_i xi*xi^T - Sum_i mi*mi^T+   = Sum_i xi*xi^T - n * mi*mi^T+   -}+   normalize (GaussianTrained distr) =+      let params (centerSum, covarianceSum, weight) =+             let c = recip weight+                 center = NC.scale c centerSum+             in  (center,+                  NC.sub (NC.scale c covarianceSum) (NC.outer center center))+      in  Gaussian $+          Array.array (fst $ Map.findMin distr, fst $ Map.findMax distr) $+          Map.toList $ fmap (gaussianParameters . params) distr++gaussian ::+   (Algo.Field prob) =>+   [(Vector prob, Matrix prob)] -> Gaussian prob+gaussian =+   consGaussian . map gaussianParameters++gaussianParameters ::+   (Algo.Field prob) =>+   (Vector prob, Matrix prob) -> (Vector prob, Matrix prob, prob)+gaussianParameters (center, covariance) =+   gaussianFromCholesky center $ Algo.chol covariance++consGaussian :: [(Vector a, Matrix a, a)] -> Gaussian a+consGaussian xs =+   Gaussian $ listArray (State 0, State $ length xs - 1) xs++gaussianFromCholesky ::+   (Algo.Field prob) =>+   Vector prob -> Matrix prob -> (Vector prob, Matrix prob, prob)+gaussianFromCholesky center covarianceChol =+   let covarianceSqrtDet = NC.prodElements $ Matrix.takeDiag covarianceChol+   in  (center, covarianceChol,+        recip (sqrt (2*pi) ^ Vector.dim center * covarianceSqrtDet))++++class CSV distr where+   toCells :: distr -> [[String]]+   parseCells :: Int -> HMMCSV.CSVParser distr++class (Ord symbol) => CSVSymbol symbol where+   cellFromSymbol :: symbol -> String+   symbolFromCell :: String -> Maybe symbol++instance CSVSymbol Char where+   cellFromSymbol = (:[])+   symbolFromCell = listToMaybe++instance CSVSymbol Int where+   cellFromSymbol = show+   symbolFromCell = maybeRead+++instance+   (Algo.Field prob, Show prob, Read prob, CSVSymbol symbol) =>+      CSV (Discrete prob symbol) where+   toCells (Discrete m) =+      map+         (\(symbol, probs) ->+            cellFromSymbol symbol : HMMCSV.cellsFromVector probs) $+      Map.toAscList m+   parseCells n =+      fmap (Discrete . Map.fromList) $+      HMMCSV.manyRowsUntilEnd $ parseSymbolProb n++parseSymbolProb ::+   (Algo.Field prob, Read prob, CSVSymbol symbol) =>+   Int -> CSV.CSVRow -> HMMCSV.CSVParser (symbol, Vector prob)+parseSymbolProb n row =+   case row of+      [] -> MT.lift $ ME.throw "missing symbol"+      c:cs ->+         liftM2 (,)+            (let str = CSV.csvFieldContent c+             in  MT.lift $ ME.fromMaybe (printf "unknown symbol %s" str) $+                 symbolFromCell str)+            (do v <- HMMCSV.parseVectorFields cs+                HMMCSV.assert (n == Vector.dim v)+                   (printf "number of states (%d) and size of probability vector (%d) mismatch"+                      n (Vector.dim v))+                return v)+++instance (Algo.Field a, Eq a, Show a, Read a) => CSV (Gaussian a) where+   toCells (Gaussian params) =+      List.intercalate [[]] $+      map+         (\(center, covarianceChol, _) ->+            HMMCSV.cellsFromVector center :+            HMMCSV.cellsFromMatrix covarianceChol) $+      Array.elems params+   parseCells n = do+      gs <- HMMCSV.manySepUntilEnd parseSingleGaussian+      HMMCSV.assert (length gs == n) $+         printf "number of states (%d) and number of Gaussians (%d) mismatch"+            n (length gs)+      return $ consGaussian gs++parseSingleGaussian ::+   (Algo.Field prob, Eq prob, Read prob) =>+   HMMCSV.CSVParser (Vector prob, Matrix prob, prob)+parseSingleGaussian = do+   center <- HMMCSV.parseNonEmptyVectorCells+   covarianceChol <- HMMCSV.parseSquareMatrixCells $ Vector.dim center+   HMMCSV.assert (isUpperTriang covarianceChol) $+      "matrices must be upper triangular"+   return $ gaussianFromCholesky center covarianceChol+++{-+Maybe this test is too strict.+It would also be ok, and certainly more intuitive+to use an orthogonal but not normalized matrix.+We could get such a matrix from the eigensystem.+-}+isUpperTriang :: (Algo.Field a, Eq a) => Matrix a -> Bool+isUpperTriang m =+   let cleared = Matrix.mapMatrixWithIndex (\(i,j) x -> if i>j then x else 0) m+   in  NC.minElement cleared == 0 &&+       NC.maxElement cleared == 0
+ src/Math/HiddenMarkovModel/Example/Circle.hs view
@@ -0,0 +1,93 @@+{- |+Example of an HMM with continuous emissions with two-dimensional observations.+We train a model to accept a parametric curve of a circle with a certain speed.+This is like "Math.HiddenMarkovModel.Example.SineWave" but in two dimensions.++The four hidden states correspond to the four quadrants.+-}+module Math.HiddenMarkovModel.Example.Circle+{-# WARNING "do not import that module, it is only intended for demonstration" #-}+   where++import qualified Math.HiddenMarkovModel as HMM+import qualified Math.HiddenMarkovModel.Distribution as Distr++import qualified Data.Packed.Matrix as Matrix+import qualified Data.Packed.Vector as Vector+import Data.Packed.Vector (Vector)++import qualified System.Random as Rnd++import qualified Control.Monad.Trans.State as MS+import Control.Monad (liftM2, replicateM)++import qualified Data.NonEmpty.Class as NonEmptyC+import qualified Data.NonEmpty as NonEmpty+import Data.Function.HT (nest)+import Data.NonEmpty ((!:))++++hmm :: HMM.Gaussian Double+hmm =+   HMM.Cons {+      HMM.initial = Vector.fromList [1/4, 1/4, 1/4, 1/4],+      HMM.transition =+         Matrix.fromLists $+            [0.9, 0.0, 0.0, 0.1] :+            [0.1, 0.9, 0.0, 0.0] :+            [0.0, 0.1, 0.9, 0.0] :+            [0.0, 0.0, 0.1, 0.9] :+            [],+      HMM.distribution =+         let cov0 = Matrix.fromLists [[0.10, -0.09],[-0.09, 0.10]]+             cov1 = Matrix.fromLists [[0.10,  0.09],[ 0.09, 0.10]]+         in  Distr.gaussian $+                (Vector.fromList [ 0.5,  0.5], cov0) :+                (Vector.fromList [-0.5,  0.5], cov1) :+                (Vector.fromList [-0.5, -0.5], cov0) :+                (Vector.fromList [ 0.5, -0.5], cov1) :+                []+   }++circleLabeled :: NonEmpty.T [] (HMM.State, Vector Double)+circleLabeled =+   NonEmpty.mapTail (take 200) $+   fmap+      (\x ->+         (HMM.state $ mod (floor (x*2/pi)) 4,+          Vector.fromList [cos x, sin x])) $+   NonEmptyC.iterate (0.1+) 0++circle :: NonEmpty.T [] (Vector Double)+circle = fmap snd circleLabeled++revealed :: NonEmpty.T [] HMM.State+revealed = HMM.reveal hmm circle++{- |+Sample multivariate normal distribution and reconstruct it from the samples.+You should obtain the same parameters.+-}+reconstructDistribution :: HMM.Gaussian Double+reconstructDistribution =+   let s0 = HMM.state 0+       gen = Distr.generate (HMM.distribution hmm) s0+   in  HMM.finishTraining $ HMM.trainSupervised 1 $ fmap ((,) s0) $+       flip MS.evalState (Rnd.mkStdGen 23) $+       liftM2 (!:) gen $ replicateM 1000 gen+++hmmTrainedSupervised :: HMM.Gaussian Double+hmmTrainedSupervised =+   HMM.finishTraining $ HMM.trainSupervised 4 circleLabeled++hmmTrainedUnsupervised :: HMM.Gaussian Double+hmmTrainedUnsupervised =+   HMM.finishTraining $ HMM.trainUnsupervised hmm circle++hmmIterativelyTrained :: HMM.Gaussian Double+hmmIterativelyTrained =+   nest 100+      (HMM.finishTraining . flip HMM.trainUnsupervised circle)+      hmm
+ src/Math/HiddenMarkovModel/Example/SineWave.hs view
@@ -0,0 +1,76 @@+{- |+Example of an HMM with continuous emissions.+We train a model to accept sine waves of a certain frequency.++There are four hidden states:+0 - rising,+1 - high,+2 - falling,+3 - low.+-}+module Math.HiddenMarkovModel.Example.SineWave+{-# WARNING "do not import that module, it is only intended for demonstration" #-}+   where++import qualified Math.HiddenMarkovModel as HMM+import qualified Math.HiddenMarkovModel.Distribution as Distr++import qualified Numeric.Container as NC+import qualified Data.Packed.Matrix as Matrix+import qualified Data.Packed.Vector as Vector++import qualified Data.NonEmpty.Class as NonEmptyC+import qualified Data.NonEmpty as NonEmpty+import Data.Function.HT (nest)+import Data.Tuple.HT (mapSnd)++++hmm :: HMM.Gaussian Double+hmm =+   HMM.Cons {+      HMM.initial = Vector.fromList [1/4, 1/4, 1/4, 1/4],+      HMM.transition =+         Matrix.fromLists $+            [0.9, 0.0, 0.0, 0.1] :+            [0.1, 0.9, 0.0, 0.0] :+            [0.0, 0.1, 0.9, 0.0] :+            [0.0, 0.0, 0.1, 0.9] :+            [],+      HMM.distribution =+         Distr.gaussian $+            (Vector.fromList [ 0], Matrix.fromLists [[1]]) :+            (Vector.fromList [ 1], Matrix.fromLists [[1]]) :+            (Vector.fromList [ 0], Matrix.fromLists [[1]]) :+            (Vector.fromList [-1], Matrix.fromLists [[1]]) :+            []+   }++sineWaveLabeled :: NonEmpty.T [] (HMM.State, Double)+sineWaveLabeled =+   NonEmpty.mapTail (take 200) $+   fmap (\x -> (HMM.state $ mod (floor (x*2/pi+0.5)) 4, sin x)) $+   NonEmptyC.iterate (0.1+) 0++sineWave :: NonEmpty.T [] Double+sineWave = fmap snd sineWaveLabeled++revealed :: NonEmpty.T [] HMM.State+revealed = HMM.reveal hmmTrainedSupervised $ fmap NC.scalar sineWave++hmmTrainedSupervised :: HMM.Gaussian Double+hmmTrainedSupervised =+   HMM.finishTraining $ HMM.trainSupervised 4 $+   fmap (mapSnd NC.scalar) sineWaveLabeled++hmmTrainedUnsupervised :: HMM.Gaussian Double+hmmTrainedUnsupervised =+   HMM.finishTraining $ HMM.trainUnsupervised hmm $ fmap NC.scalar sineWave++hmmIterativelyTrained :: HMM.Gaussian Double+hmmIterativelyTrained =+   nest 100+      (\model ->+         HMM.finishTraining $ HMM.trainUnsupervised model $+         fmap NC.scalar sineWave)+      hmm
+ src/Math/HiddenMarkovModel/Example/TrafficLight.hs view
@@ -0,0 +1,172 @@+{- |+This is an example of an HMM with discrete emissions.+We model a traffic light consisting of the colors red, yellow, green,+where only one lamp can be switched on at every point in time.+This way, when it is yellow you cannot tell immediately+whether it will switch to green or red.+We can only infer this from the light seen before.++There are four hidden states:+0 emits red, 1 emits yellow between red and green,+2 emits green, 3 emits yellow between green and red.++We quantise time in time steps.+The transition matrix of the model 'hmm' encodes+the expected duration of every state counted in time steps+and what states follow after each other.+E.g. transition probability of 0.8 of a state to itself means+that the expected duration of the state is 5 time steps (1/(1-0.8)).+However, it is a geometric distribution,+that is, shorter durations are always more probable.++The distribution of 'hmm' encodes which lights a state activates.+In our case everything is deterministic:+Every state can switch exactly one light on.++Given a sequence of observed lights+the function 'HMM.reveal' tells us the most likely sequence of states.+We test this with the light sequences in 'stateSequences'+where we already know the hidden states+as they are stored in 'labeledSequences'.+'verifyRevelation' compares the computed state sequence with the given one.++We also try some trainings in 'hmmTrainedSupervised' et.al.+-}+module Math.HiddenMarkovModel.Example.TrafficLight+{-# WARNING "do not import that module, it is only intended for demonstration" #-}+   where++import qualified Math.HiddenMarkovModel as HMM+import qualified Math.HiddenMarkovModel.Distribution as Distr++import qualified Data.Packed.Matrix as Matrix+import qualified Data.Packed.Vector as Vector++import Text.Read.HT (maybeRead)++import Control.Monad (liftM2)++import qualified Data.Map as Map+import qualified Data.NonEmpty as NonEmpty+import qualified Data.List.HT as ListHT+import Data.NonEmpty ((!:))++++data Color = Red | Yellow | Green+   deriving (Eq, Ord, Enum, Show, Read)++{- |+Using 'show' and 'read' is not always a good choice+since they must format and parse Haskell expressions+which is not of much use to the outside world.+-}+instance Distr.CSVSymbol Color where+   cellFromSymbol = show+   symbolFromCell = maybeRead+++hmm :: HMM.Discrete Double Color+hmm =+   HMM.Cons {+      HMM.initial = Vector.fromList [1/3, 1/6, 1/3, 1/6],+      HMM.transition =+         Matrix.fromLists $+            [0.8, 0.0, 0.0, 0.2] :+            [0.2, 0.8, 0.0, 0.0] :+            [0.0, 0.2, 0.8, 0.0] :+            [0.0, 0.0, 0.2, 0.8] :+            [],+      HMM.distribution =+         Distr.Discrete $ Map.fromList $+            (Red,    Vector.fromList [1,0,0,0]) :+            (Yellow, Vector.fromList [0,1,0,1]) :+            (Green,  Vector.fromList [0,0,1,0]) :+            []+   }++hmmDisturbed :: HMM.Discrete Double Color+hmmDisturbed =+   HMM.Cons {+      HMM.initial = Vector.fromList [1/4, 1/4, 1/4, 1/4],+      HMM.transition =+         Matrix.fromLists $+            [0.3, 0.2, 0.2, 0.3] :+            [0.3, 0.3, 0.2, 0.2] :+            [0.2, 0.3, 0.3, 0.2] :+            [0.2, 0.2, 0.3, 0.3] :+            [],+      HMM.distribution =+         Distr.Discrete $ Map.fromList $+            (Red,    Vector.fromList [0.6, 0.2, 0.2, 0.2]) :+            (Yellow, Vector.fromList [0.2, 0.6, 0.2, 0.6]) :+            (Green,  Vector.fromList [0.2, 0.2, 0.6, 0.2]) :+            []+   }+++red, yellowRG, green, yellowGR :: (HMM.State, Color)+red      = (HMM.state 0, Red)+yellowRG = (HMM.state 1, Yellow)+green    = (HMM.state 2, Green)+yellowGR = (HMM.state 3, Yellow)++labeledSequences :: NonEmpty.T [] (NonEmpty.T [] (HMM.State, Color))+labeledSequences =+   (red !: red : red : red :+    yellowRG : yellowRG :+    green : green : green : green : green :+    yellowGR :+    red : red : red :+    []) !:+   (green !: green : green :+    yellowGR :+    red : red : red : red :+    yellowRG :+    green : green : green : green : green :+    yellowGR : yellowGR :+    []) :+   []++{- |+Construct a Hidden Markov model by watching a set+of manually created sequences of emissions and according states.+-}+hmmTrainedSupervised :: HMM.Discrete Double Color+hmmTrainedSupervised =+   HMM.trainMany (HMM.trainSupervised 4) labeledSequences+++stateSequences :: NonEmpty.T [] (NonEmpty.T [] Color)+stateSequences = fmap (fmap snd) labeledSequences++{- |+Construct a Hidden Markov model starting from a known model+and a set of sequences that contain only the emissions, but no states.+-}+hmmTrainedUnsupervised :: HMM.Discrete Double Color+hmmTrainedUnsupervised =+   HMM.trainMany (HMM.trainUnsupervised hmm) stateSequences++{- |+Repeat unsupervised training until convergence.+-}+hmmIterativelyTrained :: HMM.Discrete Double Color+hmmIterativelyTrained =+   snd $ head $ dropWhile fst $+   ListHT.mapAdjacent (\hmm0 hmm1 -> (HMM.deviation hmm0 hmm1 > 1e-5, hmm1)) $+   iterate+      (flip HMM.trainMany stateSequences . HMM.trainUnsupervised)+      hmmDisturbed+++verifyRevelation ::+   HMM.Discrete Double Color -> NonEmpty.T [] (HMM.State, Color) -> Bool+verifyRevelation model xs =+   fmap fst xs == HMM.reveal model (fmap snd xs)++verifyRevelations :: [Bool]+verifyRevelations =+   liftM2 verifyRevelation+      [hmm, hmmDisturbed, hmmTrainedSupervised, hmmTrainedUnsupervised]+      (NonEmpty.flatten labeledSequences)
+ src/Math/HiddenMarkovModel/Named.hs view
@@ -0,0 +1,92 @@+module Math.HiddenMarkovModel.Named (+   T(..),+   Discrete,+   Gaussian,+   fromModelAndNames,+   toCSV,+   fromCSV,+   ) where++import qualified Math.HiddenMarkovModel.Distribution as Distr+import qualified Math.HiddenMarkovModel.Private as HMM+import qualified Math.HiddenMarkovModel.CSV as HMMCSV+import Math.HiddenMarkovModel.Distribution (State(..))+import Math.HiddenMarkovModel.Utility (attachOnes)++import qualified Numeric.LinearAlgebra.Algorithms as Algo+import qualified Data.Packed.Vector as Vector++import qualified Text.CSV.Lazy.String as CSV+import Text.Printf (printf)++import qualified Control.Monad.Exception.Synchronous as ME+import qualified Control.Monad.Trans.State as MS++import qualified Data.Map as Map+import qualified Data.List as List+import Data.Tuple.HT (swap)+import Data.Map (Map)+++{- |+A Hidden Markov Model with names for each state.++Although 'nameFromStateMap' and 'stateFromNameMap' are exported+you must be careful to keep them consistent when you alter them.+-}+data T distr prob =+   Cons {+      model :: HMM.T distr prob,+      nameFromStateMap :: Map State String,+      stateFromNameMap :: Map String State+   }+   deriving (Show, Read)++type Discrete prob symbol = T (Distr.Discrete prob symbol) prob+type Gaussian a = T (Distr.Gaussian a) a+++fromModelAndNames :: HMM.T distr prob -> [String] -> T distr prob+fromModelAndNames md names =+   let m = Map.fromList $ zip [State 0 ..] names+   in  Cons {+          model = md,+          nameFromStateMap = m,+          stateFromNameMap = inverseMap m+       }++inverseMap :: Map State String -> Map String State+inverseMap =+   Map.fromListWith (error "duplicate label") .+   map swap . Map.toList+++toCSV ::+   (Distr.CSV distr, Algo.Field prob, Show prob) =>+   T distr prob -> String+toCSV hmm =+   CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $+      Map.elems (nameFromStateMap hmm) : HMM.toCells (model hmm)++fromCSV ::+   (Distr.CSV distr, Algo.Field prob, Read prob) =>+   String -> ME.Exceptional String (T distr prob)+fromCSV =+   MS.evalStateT parseCSV . map HMMCSV.fixShortRow . CSV.parseCSV++parseCSV ::+   (Distr.CSV distr, Algo.Field prob, Read prob) =>+   HMMCSV.CSVParser (T distr prob)+parseCSV = do+   names <- HMMCSV.parseStringList =<< HMMCSV.getRow+   let duplicateNames =+         Map.keys $ Map.filter (> (1::Int)) $+         Map.fromListWith (+) $ attachOnes names+    in HMMCSV.assert (null duplicateNames) $+          "duplicate names: " ++ List.intercalate ", " duplicateNames+   md <- HMM.parseCSV+   let n = length names+       m = Vector.dim (HMM.initial md)+    in HMMCSV.assert (n == m) $+          printf "got %d state names for %d state" n m+   return $ fromModelAndNames md names
+ src/Math/HiddenMarkovModel/Normalized.hs view
@@ -0,0 +1,168 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{- |+Counterparts to functions in "Math.HiddenMarkovModel.Private"+that normalize interim results.+We need to do this in order to prevent+to round very small probabilities to zero.+-}+module Math.HiddenMarkovModel.Normalized where++import qualified Math.HiddenMarkovModel.Distribution as Distr+import Math.HiddenMarkovModel.Private+          (T(..), Trained(..), emission, matrixMaxMul, sumTransitions)+import Math.HiddenMarkovModel.Distribution (State(State))+import Math.HiddenMarkovModel.Utility (normalizeFactor, normalizeProb)++import qualified Numeric.Container as NC+import qualified Data.Packed.Development as Dev+import qualified Data.Packed.Vector as Vector+import Numeric.Container ((<>))+import Data.Packed.Matrix (Matrix)+import Data.Packed.Vector (Vector)++import qualified Control.Functor.HT as Functor++import qualified Data.NonEmpty.Class as NonEmptyC+import qualified Data.NonEmpty as NonEmpty+import qualified Data.Foldable as Fold+import qualified Data.List as List+import Data.Traversable (Traversable, mapAccumL)+import Data.Tuple.HT (mapFst, mapSnd, swap)+++{- |+Logarithm of the likelihood to observe the given sequence.+We return the logarithm because the likelihood can be so small+that it may be rounded to zero in the choosen number type.+-}+logLikelihood ::+   (Distr.EmissionProb distr, Floating prob,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr prob -> NonEmpty.T f emission -> prob+logLikelihood hmm = Fold.sum . fmap (log . fst) . alpha hmm++alpha ::+   (Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr prob ->+   NonEmpty.T f emission -> NonEmpty.T f (prob, Vector prob)+alpha hmm (NonEmpty.Cons x xs) =+   let normMulEmiss y = normalizeFactor . NC.mul (emission hmm y)+   in  NonEmpty.scanl+          (\(_,alphai) xi -> normMulEmiss xi (transition hmm <> alphai))+          (normMulEmiss x (initial hmm))+          xs++beta ::+   (Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f, NonEmptyC.Zip f, NonEmptyC.Reverse f) =>+   T distr prob ->+   f (prob, emission) -> NonEmpty.T f (Vector prob)+beta hmm =+   nonEmptyScanr+      (\(ci,xi) betai ->+         NC.scale (recip ci) $ NC.mul (emission hmm xi) betai <> transition hmm)+      (NC.constant 1 (NC.dim $ initial hmm))++alphaBeta ::+   (Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f, NonEmptyC.Zip f, NonEmptyC.Reverse f) =>+   T distr prob ->+   NonEmpty.T f emission ->+   (NonEmpty.T f (prob, Vector prob), NonEmpty.T f (Vector prob))+alphaBeta hmm xs =+   let calphas = alpha hmm xs+   in  (calphas,+        beta hmm $ NonEmpty.tail $ NonEmptyC.zip (fmap fst calphas) xs)+++xiFromAlphaBeta ::+   (Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f, NonEmptyC.Zip f) =>+   T distr prob ->+   NonEmpty.T f emission ->+   NonEmpty.T f (prob, Vector prob) ->+   NonEmpty.T f (Vector prob) ->+   f (Matrix prob)+xiFromAlphaBeta hmm xs calphas betas =+   let (cs,alphas) = Functor.unzip calphas+   in  NonEmptyC.zipWith4+          (\x alpha0 c1 beta1 ->+             NC.scale (recip c1) $+             NC.mul+                (NC.outer (NC.mul (emission hmm x) beta1) alpha0)+                (transition hmm))+          (NonEmpty.tail xs)+          (NonEmpty.init alphas)+          (NonEmpty.tail cs)+          (NonEmpty.tail betas)++zetaFromAlphaBeta ::+   (NC.Container Vector prob, NonEmptyC.Zip f) =>+   NonEmpty.T f (prob, Vector prob) ->+   NonEmpty.T f (Vector prob) ->+   NonEmpty.T f (Vector prob)+zetaFromAlphaBeta calphas betas =+   NonEmptyC.zipWith (NC.mul . snd) calphas betas+++{- |+Reveal the state sequence+that led most likely to the observed sequence of emissions.+It is found using the Viterbi algorithm.+-}+reveal ::+   (Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f, NonEmptyC.Reverse f) =>+   T distr prob -> NonEmpty.T f emission -> NonEmpty.T f State+reveal hmm (NonEmpty.Cons x xs) =+   fmap State $+   uncurry (nonEmptyScanr Dev.at') $+   mapFst NC.maxIndex $+   mapAccumL+      (\alphai xi ->+         swap $ mapSnd (NC.mul (emission hmm xi)) $+         matrixMaxMul (transition hmm) $ normalizeProb alphai)+      (NC.mul (emission hmm x) (initial hmm)) xs+++{- |+Variant of NonEmpty.scanr with less stack consumption.+-}+nonEmptyScanr ::+   (Traversable f, NonEmptyC.Reverse f) =>+   (a -> b -> b) -> b -> f a -> NonEmpty.T f b+nonEmptyScanr f x =+   NonEmptyC.reverse . NonEmpty.scanl (flip f) x . NonEmptyC.reverse+++{- |+Consider a superposition of all possible state sequences+weighted by the likelihood to produce the observed emission sequence.+Now train the model with respect to all of these sequences+with respect to the weights.+This is done by the Baum-Welch algorithm.+-}+trainUnsupervised ::+   (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   T distr prob -> NonEmpty.T [] emission -> Trained tdistr prob+trainUnsupervised hmm xs =+   let (alphas, betas) = alphaBeta hmm xs+       zetas = zetaFromAlphaBeta alphas betas++   in  Trained {+          trainedInitial = NonEmpty.head zetas,+          trainedTransition =+             sumTransitions hmm $ xiFromAlphaBeta hmm xs alphas betas,+          trainedDistribution =+             Distr.accumulateEmissions $ map (zip (NonEmpty.flatten xs)) $+             List.transpose $ map Vector.toList $ NonEmpty.flatten zetas+       }
+ src/Math/HiddenMarkovModel/Pattern.hs view
@@ -0,0 +1,108 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{- |+This module provides a simple way to train+the transition matrix and initial probability vector+using simple patterns of state sequences.++You may create a trained model using semigroup combinators like this:++> let a = atom $ HMM.state 0+>     b = atom $ HMM.state 1+>     distr =+>        Distr.DiscreteTrained $ Map.fromList $+>        ('a', Vector.fromList [1,2]) :+>        ('b', Vector.fromList [4,3]) :+>        ('c', Vector.fromList [0,1]) :+>        []+> in  finish 2 distr $ replicate 5 $ replicate 10 a <> replicate 20 b+-}+module Math.HiddenMarkovModel.Pattern (+   T,+   atom,+   append,+   replicate,+   finish,+   ) where++import qualified Math.HiddenMarkovModel.Distribution as Distr+import qualified Math.HiddenMarkovModel as HMM+import Math.HiddenMarkovModel.Private (Trained(..))+import Math.HiddenMarkovModel.Distribution (State(State))++import qualified Numeric.LinearAlgebra.Algorithms as Algo+import qualified Numeric.Container as NC+import qualified Data.Packed.Vector as Vector+import Data.Packed.Matrix (Matrix)+import Data.Packed.Vector (Vector)++import qualified Data.Map as Map+import Data.Semigroup (Semigroup, (<>), times1p)++import Prelude hiding (replicate)+++newtype T prob = Cons (Int -> (State, Matrix prob, State))++atom ::+   (NC.Container Vector prob) =>+   State -> T prob+atom s = Cons $ \n -> (s, NC.konst 0 (n,n), s)+++instance (Algo.Field prob) => Semigroup (T prob) where+   (<>) = append+   times1p k = replicate $ fromIntegral (k-1)+++infixl 5 `append`++append ::+   (NC.Container Vector prob) =>+   T prob -> T prob -> T prob+append (Cons f) (Cons g) =+   Cons $ \n ->+      case (f n, g n) of+         ((sai, ma, sao), (sbi, mb, sbo)) ->+            (sai, increment (sbi,sao) 1 $ NC.add ma mb, sbo)++replicate ::+   (NC.Container Vector prob) =>+   Int -> T prob -> T prob+replicate ki (Cons f) =+   Cons $ \n ->+      case f n of+         (si, m, so) ->+            let k = fromIntegral ki+            in  (si, increment (si,so) (k-1) $ NC.scale k m, so)++increment ::+   (NC.Container Vector a) =>+   (State, State) -> a -> Matrix a -> Matrix a+increment (State i, State j) x m  =  NC.accum m (+) [((i,j), x)]+++finish ::+   (NC.Container Vector prob) =>+   Int -> tdistr -> T prob -> Trained tdistr prob+finish n tdistr (Cons f) =+   case f n of+      (State si, m, _so) ->+         Trained {+            trainedInitial = NC.assoc n 0 [(si,1)],+            trainedTransition = m,+            trainedDistribution = tdistr+         }+++_example :: HMM.DiscreteTrained Double Char+_example =+   let a = atom $ HMM.state 0+       b = atom $ HMM.state 1+       distr =+          Distr.DiscreteTrained $ Map.fromList $+          ('a', Vector.fromList [1,2]) :+          ('b', Vector.fromList [4,3]) :+          ('c', Vector.fromList [0,1]) :+          []+   in  finish 2 distr $ replicate 5 $ replicate 10 a <> replicate 20 b
+ src/Math/HiddenMarkovModel/Private.hs view
@@ -0,0 +1,286 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+module Math.HiddenMarkovModel.Private where++import qualified Math.HiddenMarkovModel.Distribution as Distr+import qualified Math.HiddenMarkovModel.CSV as HMMCSV+import Math.HiddenMarkovModel.Distribution (State(State))++import qualified Numeric.LinearAlgebra.Algorithms as Algo+import qualified Numeric.LinearAlgebra.Util as LinAlg+import qualified Numeric.Container as NC+import qualified Data.Packed.Development as Dev+import qualified Data.Packed.Matrix as Matrix+import qualified Data.Packed.Vector as Vector+import Numeric.Container ((<>))+import Data.Packed.Matrix (Matrix)+import Data.Packed.Vector (Vector)++import qualified Data.NonEmpty.Class as NonEmptyC+import qualified Data.NonEmpty as NonEmpty+import qualified Data.Semigroup as Sg+import qualified Data.List as List+import Data.Traversable (Traversable, mapAccumL)+import Data.Tuple.HT (mapPair, mapFst, mapSnd, swap)+++{- |+A Hidden Markov model consists of a number of (hidden) states+and a set of emissions.+There is a vector for the initial probability of each state+and a matrix containing the probability for switching+from one state to another one.+The 'distribution' field points to probability distributions+that associate every state with emissions of different probability.+Famous distribution instances are discrete and Gaussian distributions.+See "Math.HiddenMarkovModel.Distribution" for details.++The transition matrix is transposed+with respect to popular HMM descriptions.+But I think this is the natural orientation, because this way+you can write \"transition matrix times probability column vector\".++The type has two type parameters,+although the one for the distribution would be enough.+However, replacing @prob@ by @Distr.Probability distr@+would prohibit the derived Show and Read instances.+-}+data T distr prob =+   Cons {+      initial :: Vector prob,+      transition :: Matrix prob,+      distribution :: distr+   }+   deriving (Show, Read)+++emission ::+   (Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   T distr prob -> emission -> Vector prob+emission  =  Distr.emissionProb . distribution+++forward ::+   (Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr prob -> NonEmpty.T f emission -> prob+forward hmm = NC.sumElements . NonEmpty.last . alpha hmm++alpha ::+   (Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr prob ->+   NonEmpty.T f emission -> NonEmpty.T f (Vector prob)+alpha hmm (NonEmpty.Cons x xs) =+   NonEmpty.scanl+      (\alphai xi -> NC.mul (emission hmm xi) (transition hmm <> alphai))+      (NC.mul (emission hmm x) (initial hmm))+      xs+++backward ::+   (Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr prob -> NonEmpty.T f emission -> prob+backward hmm (NonEmpty.Cons x xs) =+   NC.sumElements $+   NC.mul (initial hmm) $+   NC.mul (emission hmm x) $+   NonEmpty.head $ beta hmm xs++beta ::+   (Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr prob ->+   f emission -> NonEmpty.T f (Vector prob)+beta hmm =+   NonEmpty.scanr+      (\xi betai -> NC.mul (emission hmm xi) betai <> transition hmm)+      (NC.constant 1 (NC.dim $ initial hmm))+++alphaBeta ::+   (Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr prob ->+   NonEmpty.T f emission ->+   (prob, NonEmpty.T f (Vector prob), NonEmpty.T f (Vector prob))+alphaBeta hmm xs =+   let alphas = alpha hmm xs+       betas = beta hmm $ NonEmpty.tail xs+       recipLikelihood = recip $ NC.sumElements $ NonEmpty.last alphas+   in  (recipLikelihood, alphas, betas)++++xiFromAlphaBeta ::+   (Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   T distr prob -> prob ->+   NonEmpty.T [] emission ->+   NonEmpty.T [] (Vector prob) ->+   NonEmpty.T [] (Vector prob) ->+   [Matrix prob]+xiFromAlphaBeta hmm recipLikelihood xs alphas betas =+   zipWith3+      (\x alpha0 beta1 ->+         NC.scale recipLikelihood $+         NC.mul+            (NC.outer (NC.mul (emission hmm x) beta1) alpha0)+            (transition hmm))+      (NonEmpty.tail xs)+      (NonEmpty.init alphas)+      (NonEmpty.tail betas)++zetaFromXi ::+   (Distr.Probability distr ~ prob, Num prob, NC.Product prob) =>+   T distr prob -> [Matrix prob] -> [Vector prob]+zetaFromXi hmm xis =+   map (NC.constant 1 (Matrix.rows $ transition hmm) <>) xis++zetaFromAlphaBeta ::+   (NC.Container Vector prob) =>+   prob ->+   NonEmpty.T [] (Vector prob) ->+   NonEmpty.T [] (Vector prob) ->+   NonEmpty.T [] (Vector prob)+zetaFromAlphaBeta recipLikelihood alphas betas =+   fmap (NC.scale recipLikelihood) $+   NonEmptyC.zipWith NC.mul alphas betas+++{- |+In constrast to Math.HiddenMarkovModel.reveal+this does not normalize the vector.+This is slightly simpler but for long sequences+the product of probabilities might be smaller+than the smallest representable number.+-}+reveal ::+   (Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr prob -> NonEmpty.T f emission -> NonEmpty.T f State+reveal hmm (NonEmpty.Cons x xs) =+   fmap State $+   uncurry (NonEmpty.scanr Dev.at') $+   mapFst NC.maxIndex $+   mapAccumL+      (\alphai xi ->+         swap $ mapSnd (NC.mul (emission hmm xi)) $+         matrixMaxMul (transition hmm) alphai)+      (NC.mul (emission hmm x) (initial hmm)) xs++matrixMaxMul ::+   (NC.Container Vector a) =>+   Matrix a -> Vector a -> (Vector Int, Vector a)+matrixMaxMul m v =+   mapPair (Vector.fromList, Vector.fromList) $ unzip $+   map ((\x -> (NC.maxIndex x, NC.maxElement x)) . NC.mul v) $+   Matrix.toRows m++++{- |+A trained model is a temporary form of a Hidden Markov model+that we need during the training on multiple training sequences.+It allows to collect knowledge over many sequences with 'mergeTrained',+even with mixed supervised and unsupervised training.+You finish the training by converting the trained model+back to a plain modul using 'finishTraining'.++You can create a trained model in three ways:++* supervised training using an emission sequence with associated states,++* unsupervised training using an emission sequence and an existing Hidden Markov Model,++* derive it from state sequence patterns, cf. "Math.HiddenMarkovModel.Pattern".+-}+data Trained distr prob =+   Trained {+      trainedInitial :: Vector prob,+      trainedTransition :: Matrix prob,+      trainedDistribution :: distr+   }+   deriving (Show, Read)+++sumTransitions ::+   (NC.Container Vector e, Num e) =>+   T distr t -> [Matrix e] -> Matrix e+sumTransitions hmm =+   foldl NC.add (NC.konst 0 $ LinAlg.size $ transition hmm)++{- |+Baum-Welch algorithm+-}+trainUnsupervised ::+   (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   T distr prob -> NonEmpty.T [] emission -> Trained tdistr prob+trainUnsupervised hmm xs =+   let (recipLikelihood, alphas, betas) = alphaBeta hmm xs+       zetas = zetaFromAlphaBeta recipLikelihood alphas betas++   in  Trained {+          trainedInitial = NonEmpty.head zetas,+          trainedTransition =+             sumTransitions hmm $+             xiFromAlphaBeta hmm recipLikelihood xs alphas betas,+          trainedDistribution =+             Distr.accumulateEmissions $ map (zip (NonEmpty.flatten xs)) $+             List.transpose $ map Vector.toList $ NonEmpty.flatten zetas+       }+++mergeTrained ::+   (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   Trained tdistr prob -> Trained tdistr prob -> Trained tdistr prob+mergeTrained hmm0 hmm1 =+   Trained {+      trainedInitial = NC.add (trainedInitial hmm0) (trainedInitial hmm1),+      trainedTransition =+         NC.add (trainedTransition hmm0) (trainedTransition hmm1),+      trainedDistribution =+         Distr.combine+            (trainedDistribution hmm0) (trainedDistribution hmm1)+   }++instance+   (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,+    Distr.Probability distr ~ prob) =>+      Sg.Semigroup (Trained tdistr prob) where+   (<>) = mergeTrained+++toCells ::+   (Distr.CSV distr, Algo.Field prob, Show prob) =>+   T distr prob -> [[String]]+toCells hmm =+   (HMMCSV.cellsFromVector $ initial hmm) :+   (HMMCSV.cellsFromMatrix $ transition hmm) +++   [] :+   (Distr.toCells $ distribution hmm)++parseCSV ::+   (Distr.CSV distr, Algo.Field prob, Read prob) =>+   HMMCSV.CSVParser (T distr prob)+parseCSV = do+   v <- HMMCSV.parseNonEmptyVectorCells+   m <- HMMCSV.parseSquareMatrixCells $ Vector.dim v+   HMMCSV.skipEmptyRow+   distr <- Distr.parseCells $ Vector.dim v+   return $ Cons {+      initial = v,+      transition = m,+      distribution = distr+   }
+ src/Math/HiddenMarkovModel/Test.hs view
@@ -0,0 +1,143 @@+module Math.HiddenMarkovModel.Test where++import qualified Math.HiddenMarkovModel as HMM+import qualified Math.HiddenMarkovModel.Normalized as Normalized+import qualified Math.HiddenMarkovModel.Private as Priv+import qualified Math.HiddenMarkovModel.Distribution as Distr++import qualified Numeric.Container as NC+import qualified Data.Packed.Matrix as Matrix+import qualified Data.Packed.Vector as Vector+import Data.Packed.Matrix (Matrix)+import Data.Packed.Vector (Vector)++import qualified System.Random as Rnd++import qualified Data.NonEmpty.Class as NonEmptyC+import qualified Data.NonEmpty as NonEmpty+import qualified Data.Foldable as Fold+import qualified Data.Map as Map+import Data.NonEmpty ((!:))+++hmm :: HMM.Discrete Double Char+hmm =+   HMM.Cons {+      HMM.initial = Vector.fromList [0.1, 0.2, 0.3, 0.4],+      HMM.transition =+         Matrix.fromLists $+            [0.7, 0.1, 0.0, 0.2] :+            [0.1, 0.6, 0.1, 0.0] :+            [0.1, 0.2, 0.7, 0.0] :+            [0.1, 0.1, 0.2, 0.8] :+            [],+      HMM.distribution =+         Distr.Discrete $ Map.fromList $+            ('a', Vector.fromList [1,0,0,0]) :+            ('b', Vector.fromList [0,1,0,1]) :+            ('c', Vector.fromList [0,0,1,0]) :+            []+   }+++sequ :: NonEmpty.T [] Char+sequ = 'a' !: take 20 (HMM.generate hmm (Rnd.mkStdGen 42))++{- |+Should all be equal.+-}+sequLikelihood :: ((Double, Double), Double, NonEmpty.T [] Double)+sequLikelihood =+   ((Priv.forward hmm sequ, Priv.backward hmm sequ),+    exp $ Normalized.logLikelihood hmm sequ,+    NonEmptyC.zipWith NC.dot+       (Priv.alpha hmm sequ)+       (Priv.beta hmm $ NonEmpty.tail sequ))++{- |+Should all be one.+-}+sequLikelihoodNormalized :: NonEmpty.T [] Double+sequLikelihoodNormalized =+   let (calphas,betas) = Normalized.alphaBeta hmm sequ+   in  NonEmptyC.zipWith NC.dot (fmap snd calphas) betas+++{- |+Lists should be equal, but the first list contains one less element.+-}+zetas ::+   ([Vector Double],+    NonEmpty.T [] (Vector Double),+    NonEmpty.T [] (Vector Double))+zetas =+   let (recipLikelihood, alphas, betas) = Priv.alphaBeta hmm sequ+   in  (Priv.zetaFromXi hmm $+           Priv.xiFromAlphaBeta hmm recipLikelihood sequ alphas betas,+        Priv.zetaFromAlphaBeta recipLikelihood alphas betas,+        uncurry Normalized.zetaFromAlphaBeta $+        Normalized.alphaBeta hmm sequ)++{- |+Quick test of zetas - result should be @(True, very small, very small)@.+-}+zetasDiff :: (Bool, Double, Double)+zetasDiff =+   case zetas of+      (z0,z1,z2) ->+         (length z0 == length (NonEmpty.tail z1) &&+          length z0 == length (NonEmpty.tail z2),+          maximum $ map NC.normInf $ zipWith NC.sub z0 $ NonEmpty.init z1,+          NonEmpty.maximum $ fmap NC.normInf $ NonEmptyC.zipWith NC.sub z1 z2)++{- |+Lists should be equal+-}+xis :: ([Matrix Double], [Matrix Double])+xis =+   let (recipLikelihood, alphas, betas) = Priv.alphaBeta hmm sequ+   in  (Priv.xiFromAlphaBeta hmm recipLikelihood sequ alphas betas,+        uncurry (Normalized.xiFromAlphaBeta hmm sequ) $+        Normalized.alphaBeta hmm sequ)++{- |+Quick test of xis - result should be @(True, very small)@.+-}+xisDiff :: (Bool, Double)+xisDiff =+   case xis of+      (x0,x1) ->+         (length x0 == length x1,+          maximum $ map (NC.normInf . Matrix.flatten) $ zipWith NC.sub x0 x1)+++reveal :: Bool+reveal =+   Normalized.reveal hmm sequ == Priv.reveal hmm sequ+++trainUnsupervised ::+   (HMM.DiscreteTrained Double Char,+    HMM.DiscreteTrained Double Char)+trainUnsupervised =+   (Priv.trainUnsupervised hmm sequ,+    Normalized.trainUnsupervised hmm sequ)++trainUnsupervisedDiff :: (Double, Double, (Bool, Double))+trainUnsupervisedDiff =+   case trainUnsupervised of+      (hmm0,hmm1) ->+         (NC.normInf $ Matrix.flatten $ NC.sub+             (Priv.trainedTransition hmm0) (Priv.trainedTransition hmm1),+          NC.normInf $ NC.sub+             (Priv.trainedInitial hmm0) (Priv.trainedInitial hmm1),+          case (Priv.trainedDistribution hmm0, Priv.trainedDistribution hmm1) of+             (Distr.DiscreteTrained m0, Distr.DiscreteTrained m1) ->+                (Map.size m0 == Map.size m1,+                 Fold.maximum $ fmap NC.normInf $+                    Map.intersectionWith NC.sub m0 m1))+++nonEmptyScanr :: Int -> [Int] -> Bool+nonEmptyScanr x xs =+   Normalized.nonEmptyScanr (-) x xs == NonEmpty.scanr (-) x xs
+ src/Math/HiddenMarkovModel/Utility.hs view
@@ -0,0 +1,35 @@+{-# LANGUAGE FlexibleContexts #-}+module Math.HiddenMarkovModel.Utility where++import qualified Numeric.Container as NC+import Data.Packed.Vector (Vector)++import qualified System.Random as Rnd++import qualified Control.Monad.Trans.State as MS+++normalizeProb ::+   (NC.Container Vector a, Fractional a) => Vector a -> Vector a+normalizeProb = snd . normalizeFactor++normalizeFactor ::+   (NC.Container Vector a, Fractional a) =>+   Vector a -> (a, Vector a)+normalizeFactor xs =+   let c = NC.sumElements xs+   in  (c, NC.scale (recip c) xs)++-- see htam:Stochastic+randomItemProp ::+   (Rnd.RandomGen g, Rnd.Random b, Num b, Ord b) =>+   [(a,b)] -> MS.State g a+randomItemProp props =+   let (keys,ps) = unzip props+   in  do p <- MS.state (Rnd.randomR (0, sum ps))+          return $+             fst $ head $ dropWhile ((0<=) . snd) $+             zip keys $ tail $ scanl (-) p ps++attachOnes :: (Num b) => [a] -> [(a,b)]+attachOnes = map (flip (,) 1)