hmm-lapack (empty) → 0.3

raw patch · 19 files changed

+2526/−0 lines, 19 filesdep +QuickCheckdep +basedep +boxessetup-changed

Dependencies added: QuickCheck, base, boxes, comfort-array, containers, deepseq, explicit-exception, fixed-length, hmm-lapack, lapack, lazy-csv, netlib-ffi, non-empty, prelude-compat, random, semigroups, tfp, transformers, utility-ht

Files

• Changes.md +3/−0
• LICENSE +30/−0
• Setup.lhs +3/−0
• hmm-lapack.cabal +101/−0
• src/Math/HiddenMarkovModel.hs +221/−0
• src/Math/HiddenMarkovModel/CSV.hs +160/−0
• src/Math/HiddenMarkovModel/Distribution.hs +481/−0
• src/Math/HiddenMarkovModel/Example/Circle.hs +12/−0
• src/Math/HiddenMarkovModel/Example/CirclePrivate.hs +123/−0
• src/Math/HiddenMarkovModel/Example/SineWave.hs +90/−0
• src/Math/HiddenMarkovModel/Example/TrafficLight.hs +50/−0
• src/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs +164/−0
• src/Math/HiddenMarkovModel/Named.hs +117/−0
• src/Math/HiddenMarkovModel/Normalized.hs +175/−0
• src/Math/HiddenMarkovModel/Pattern.hs +121/−0
• src/Math/HiddenMarkovModel/Private.hs +335/−0
• src/Math/HiddenMarkovModel/Test.hs +258/−0
• src/Math/HiddenMarkovModel/Utility.hs +72/−0
• test/Main.hs +10/−0

@@ -0,0 +1,3 @@+## 0.1++* `Distribution.Estimate` turned into a multi-parameter type class.

@@ -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.

@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain

@@ -0,0 +1,101 @@+Name:                hmm-lapack+Version:             0.3+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://hub.darcs.net/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+Extra-Source-Files:+  Changes.md++Source-Repository this+  Tag:         0.3+  Type:        darcs+  Location:    http://hub.darcs.net/thielema/hmm-hmatrix++Source-Repository head+  Type:        darcs+  Location:    http://hub.darcs.net/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+    Math.HiddenMarkovModel.Test+  Other-Modules:+    Math.HiddenMarkovModel.Example.TrafficLightPrivate+    Math.HiddenMarkovModel.Example.CirclePrivate+    Math.HiddenMarkovModel.Normalized+    Math.HiddenMarkovModel.Private+    Math.HiddenMarkovModel.Utility+    Math.HiddenMarkovModel.CSV+  Build-Depends:+    lapack >=0.2 && <0.3,+    fixed-length >=0.2 && <0.3,+    tfp >=1.0 && <1.1,+    netlib-ffi >=0.1.1 && <0.2,+    comfort-array >=0.2 && <0.3,+    QuickCheck >=2.5 && <3,+    explicit-exception >=0.1.7 && <0.2,+    boxes >=0.1.5 && <0.2,+    lazy-csv >=0.5 && <0.6,+    random >=1.0 && <1.2,+    transformers >= 0.2 && <0.6,+    non-empty >=0.2.1 && <0.4,+    semigroups >=0.17 && <0.19,+    containers >=0.4.2 && <0.7,+    utility-ht >=0.0.12 && <0.1,+    deepseq >=1.3 && <1.5,+    prelude-compat >=0.0 && <0.1,+    base >=4.5 && <5+  Hs-Source-Dirs:      src+  Default-Language:    Haskell2010+  GHC-Options:         -Wall++Test-Suite hmm-test+  Type: exitcode-stdio-1.0+  Build-Depends:+    hmm-lapack,+    QuickCheck,+    base+  Main-Is: Main.hs+  Hs-Source-Dirs:      test+  Default-Language:    Haskell2010+  GHC-Options:         -Wall

@@ -0,0 +1,221 @@+{-# LANGUAGE TypeFamilies #-}+module Math.HiddenMarkovModel (+   T(..),+   Discrete, DiscreteTrained,+   Gaussian, GaussianTrained,+   uniform,+   generate,+   generateLabeled,+   probabilitySequence,+   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.Utility+          (SquareMatrix, squareConstant,+           randomItemProp, normalizeProb, attachOnes)++import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import qualified Numeric.LAPACK.Scalar as Scalar++import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable as StorableArray+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.Array.Comfort.Boxed as Array++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 Data.Traversable (Traversable, mapAccumL)+import Data.Foldable (Foldable)++++type DiscreteTrained symbol sh prob =+         Trained (Distr.DiscreteTrained symbol sh prob) sh prob+type Discrete symbol sh prob = T (Distr.Discrete symbol sh prob) sh prob++type GaussianTrained emiSh stateSh a =+         Trained (Distr.GaussianTrained emiSh stateSh a) stateSh a+type Gaussian emiSh stateSh a = T (Distr.Gaussian emiSh stateSh a) stateSh 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.StateShape distr ~ sh, Shape.C sh,+    Distr.Probability distr ~ prob) =>+   distr -> T distr sh prob+uniform distr =+   let sh = Distr.statesShape distr+       c = recip $ fromIntegral $ Shape.size sh+   in  Cons {+          initial = Vector.constant sh c,+          transition = squareConstant sh c,+          distribution = distr+       }+++probabilitySequence ::+   (Traversable f, Distr.EmissionProb distr,+    Distr.StateShape distr ~ sh, Shape.Indexed sh, Shape.Index sh ~ state,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   T distr sh prob -> f (state, emission) -> f prob+probabilitySequence hmm =+   snd+   .+   mapAccumL+      (\index (s, e) ->+         ((transition hmm StorableArray.!) . flip (,) s,+          index s * Distr.emissionStateProb (distribution hmm) e s))+      (initial hmm StorableArray.!)++generate ::+   (Rnd.RandomGen g, Ord prob, Rnd.Random prob, Distr.Generate distr,+    Distr.StateShape distr ~ sh, Shape.Indexed sh, Shape.Index sh ~ state,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   T distr sh prob -> g -> [emission]+generate hmm = map snd . generateLabeled hmm++generateLabeled ::+   (Rnd.RandomGen g, Ord prob, Rnd.Random prob, Distr.Generate distr,+    Distr.StateShape distr ~ sh, Shape.Indexed sh, Shape.Index sh ~ state,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   T distr sh prob -> g -> [(state, emission)]+generateLabeled hmm =+   MS.evalState $+   flip MS.evalStateT (initial hmm) $+   Monad.repeat $ MS.StateT $ \v0 -> do+      s <-+         randomItemProp $+         zip (Shape.indices $ StorableArray.shape v0) (Vector.toList v0)+      x <- Distr.generate (distribution hmm) s+      return ((s, x), Matrix.takeColumn (transition hmm) s)++++{- |+Contribute a manually labeled emission sequence to a HMM training.+-}+trainSupervised ::+   (Distr.StateShape distr ~ sh, Shape.Index sh ~ state,+    Distr.Estimate tdistr distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   sh -> NonEmpty.T [] (state, emission) -> Trained tdistr sh prob+trainSupervised sh xs =+   let getState (s, _x) = s+   in  Trained {+          trainedInitial =+             StorableArray.fromAssociations sh 0+                [(getState (NonEmpty.head xs), 1)],+          trainedTransition =+             Matrix.transpose $+             StorableArray.accumulate (+) (squareConstant sh 0) $+             attachOnes $ NonEmpty.mapAdjacent (,) $ fmap getState xs,+          trainedDistribution =+             Distr.accumulateEmissions $ Array.map attachOnes $+             Array.accumulate (flip (:))+                (Array.fromList sh $ replicate (Shape.size sh) [])+                (NonEmpty.flatten xs)+       }++finishTraining ::+   (Shape.C sh, Eq sh,+    Distr.Estimate tdistr distr, Distr.Probability distr ~ prob) =>+   Trained tdistr sh prob -> T distr sh prob+finishTraining hmm =+   Cons {+      initial = normalizeProb $ trainedInitial hmm,+      transition = normalizeProbColumns $ trainedTransition hmm,+      distribution = Distr.normalize $ trainedDistribution hmm+   }++normalizeProbColumns ::+   (Shape.C sh, Eq sh, Class.Real a) => SquareMatrix sh a -> SquareMatrix sh a+normalizeProbColumns m =+   Matrix.scaleColumns (StorableArray.map recip (Matrix.columnSums m)) m++trainMany ::+   (Shape.C sh, Eq sh,+    Distr.Estimate tdistr distr, Distr.Probability distr ~ prob,+    Foldable f) =>+   (trainingData -> Trained tdistr sh prob) ->+   NonEmpty.T f trainingData -> T distr sh 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 ::+   (Shape.InvIndexed sh, Eq sh, Class.Real prob, Ord prob) =>+   T distr sh prob -> T distr sh prob -> prob+deviation hmm0 hmm1 =+   deviationVec (initial hmm0) (initial hmm1)+   `max`+   deviationVec (transition hmm0) (transition hmm1)++deviationVec ::+   (Shape.InvIndexed sh, Eq sh, Class.Real a) =>+   StorableArray.Array sh a -> StorableArray.Array sh a -> a+deviationVec =+   getDeviation $ Class.switchReal deviationVecAux deviationVecAux++newtype Deviation f a = Deviation {getDeviation :: f a -> f a -> a}++deviationVecAux ::+   (Shape.InvIndexed sh, Eq sh, Ord a, Class.Real a, Scalar.RealOf a ~ a) =>+   Deviation (StorableArray.Array sh) a+deviationVecAux =+   Deviation $ \x y ->+      Scalar.absolute $ snd $ Vector.argAbsMaximum $ Vector.sub x y+++toCSV ::+   (Distr.ToCSV distr, Shape.Indexed sh, Class.Real prob, Show prob) =>+   T distr sh prob -> String+toCSV hmm =+   CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $+   toCells hmm++fromCSV ::+   (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh,+    Shape.Indexed stateSh, Shape.Index stateSh ~ state,+    Class.Real prob, Read prob) =>+   (Int -> stateSh) -> String -> ME.Exceptional String (T distr stateSh prob)+fromCSV makeShape =+   MS.evalStateT (parseCSV makeShape) . map HMMCSV.fixShortRow . CSV.parseCSV

@@ -0,0 +1,160 @@+module Math.HiddenMarkovModel.CSV where++import Math.HiddenMarkovModel.Utility (SquareMatrix, vectorDim)++import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix (ZeroInt)+import Numeric.LAPACK.Vector (Vector)++import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable as ComfortArray+import qualified Data.Array.Comfort.Shape as Shape++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.Reverse.StrictElement as Rev+import qualified Data.List.HT as ListHT+++cellsFromVector ::+   (Shape.C sh, Show a, Class.Real a) => Vector sh a -> [String]+cellsFromVector = map show . Vector.toList++cellsFromSquare ::+   (Shape.Indexed sh, Show a, Class.Real a) => SquareMatrix sh a -> [[String]]+cellsFromSquare = map (map show . Vector.toList) . Matrix.toRows++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, Class.Real a) =>+   CSVParser (Vector ZeroInt a)+parseVectorCells =+   parseVectorFields =<< getRow++-- ToDo: Maybe check row consistency already here?+parseVectorFields ::+   (Read a, Class.Real a) =>+   CSV.CSVRow -> CSVParser (Vector ZeroInt a)+parseVectorFields =+   MT.lift . fmap Vector.autoFromList . mapM parseNumberCell .+   Rev.dropWhile (null . CSV.csvFieldContent)++parseNonEmptyVectorCells ::+   (Read a, Class.Real a) =>+   CSVParser (Vector ZeroInt a)+parseNonEmptyVectorCells = do+   v <- parseVectorCells+   assert (vectorDim 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 ::+   (Shape.C sh, Read a, Class.Real a) =>+   sh -> CSVParser (SquareMatrix sh a)+parseSquareMatrixCells sh = do+   let n = Shape.size sh+   rows <- replicateM n parseVectorCells+   assert (not $ null rows) "no rows"+   assert (all ((n==) . vectorDim) rows) "inconsistent matrix dimensions"+   return $+      ComfortArray.reshape (MatrixShape.square MatrixShape.RowMajor sh) $+      Matrix.fromRows (Shape.ZeroBased n) rows++parseStringList :: CSV.CSVRow -> CSVParser [String]+parseStringList =+   MT.lift . mapM cellContent .+   Rev.dropWhile (null . CSV.csvFieldContent)

@@ -0,0 +1,481 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE MultiParamTypeClasses #-}+module Math.HiddenMarkovModel.Distribution (+   Emission, Probability, StateShape,+   Info(..), Generate(..), EmissionProb(..), Estimate(..),++   Discrete(..), DiscreteTrained(..),+   Gaussian(..), GaussianTrained(..), gaussian,++   ToCSV(..), FromCSV(..), HMMCSV.CSVParser, CSVSymbol(..),+   ) where++import qualified Math.HiddenMarkovModel.CSV as HMMCSV+import Math.HiddenMarkovModel.Utility (SquareMatrix, randomItemProp, vectorDim)++import qualified Numeric.LAPACK.Matrix.HermitianPositiveDefinite as HermitianPD+import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian+import qualified Numeric.LAPACK.Matrix.Triangular as Triangular+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix ((<#))+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Format (FormatArray, Format(format))++import qualified Numeric.Netlib.Class as Class+import Foreign.Storable (Storable)++import qualified Data.Array.Comfort.Storable as StorableArray+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.Array.Comfort.Boxed as Array+import Data.Array.Comfort.Boxed (Array, (!))++import qualified System.Random as Rnd++import qualified Text.CSV.Lazy.String as CSV+import qualified Text.PrettyPrint.Boxes as TextBox+import Text.PrettyPrint.Boxes ((<>), (<+>))+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.DeepSeq (NFData, rnf)+import Control.Monad (liftM2)+import Control.Applicative (liftA2, (<|>))++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.List.HT as ListHT+import qualified Data.List as List+import Data.Functor.Identity (Identity(Identity), runIdentity)+import Data.Foldable (Foldable, foldMap)+import Data.Tuple.HT (mapFst, fst3, swap)+import Data.Monoid (Endo(Endo, appEndo))+import Data.Map (Map)+import Data.Maybe (fromMaybe, listToMaybe)++import Prelude ()+import Prelude2010+++type HermitianMatrix sh = Hermitian.Hermitian sh+type UpperTriangular sh = Triangular.Upper sh+++type family Probability distr+type family Emission distr+type family StateShape distr+++class (Class.Real (Probability distr)) => Info distr where+   statesShape :: distr -> StateShape distr++class (Class.Real (Probability distr)) => Generate distr where+   generate ::+      (Rnd.RandomGen g, Emission distr ~ emission, StateShape distr ~ sh) =>+      distr -> Shape.Index sh -> MS.State g emission++class+   (Shape.Indexed (StateShape distr), Class.Real (Probability distr)) =>+      EmissionProb distr where+   {-+   This function could be implemented generically in terms of emissionStateProb+   but that would require an Info constraint.+   -}+   emissionProb ::+      distr -> Emission distr -> Vector (StateShape distr) (Probability distr)+   emissionStateProb ::+      distr -> Emission distr -> Shape.Index (StateShape distr) -> Probability distr+   emissionStateProb distr e s = emissionProb distr e StorableArray.! s++class+   (Distribution tdistr ~ distr, Trained distr ~ tdistr, EmissionProb distr) =>+      Estimate tdistr distr where+   type Distribution tdistr+   type Trained distr+   accumulateEmissions ::+      (Probability distr ~ prob, StateShape distr ~ sh) =>+      Array sh [(Emission distr, prob)] -> tdistr+   -- could as well be in Semigroup class+   combine :: tdistr -> tdistr -> tdistr+   normalize :: tdistr -> distr+++newtype Discrete symbol sh prob = Discrete (Map symbol (Vector sh prob))+   deriving (Show)++newtype+   DiscreteTrained symbol sh prob =+      DiscreteTrained (Map symbol (Vector sh prob))+   deriving (Show)++type instance Probability (Discrete symbol sh prob) = prob+type instance Emission (Discrete symbol sh prob) = symbol+type instance StateShape (Discrete symbol sh prob) = sh+++instance+   (NFData sh, NFData prob, NFData symbol) =>+      NFData (Discrete symbol sh prob) where+   rnf (Discrete m) = rnf m++instance+   (NFData sh, NFData prob, NFData symbol) =>+      NFData (DiscreteTrained symbol sh prob) where+   rnf (DiscreteTrained m) = rnf m++instance+   (FormatArray sh, Class.Real prob, Format symbol) =>+      Format (Discrete symbol sh prob) where+   format fmt (Discrete m) =+      TextBox.vsep 1 TextBox.left $+      map (\(sym,v) -> format fmt sym <> TextBox.char ':' <+> format fmt v) $+      Map.toAscList m++instance+   (Shape.C sh, Class.Real prob, Ord symbol) =>+      Info (Discrete symbol sh prob) where+   statesShape (Discrete m) = StorableArray.shape $ snd $ Map.findMin m++instance+   (Shape.Indexed sh, Class.Real prob, Ord symbol, Ord prob, Rnd.Random prob) =>+      Generate (Discrete symbol sh prob) where+   generate (Discrete m) state =+      randomItemProp $ Map.toAscList $ fmap (StorableArray.! state) m++instance+   (Shape.Indexed sh, Class.Real prob, Ord symbol) =>+      EmissionProb (Discrete symbol sh prob) where+   emissionProb (Discrete m) =+      mapLookup "emitDiscrete: unknown emission symbol" m++instance+   (Shape.Indexed sh, Eq sh, Class.Real prob, Ord symbol) =>+      Estimate (DiscreteTrained symbol sh prob) (Discrete symbol sh prob) where+   type Distribution (DiscreteTrained symbol sh prob) = Discrete symbol sh prob+   type Trained (Discrete symbol sh prob) = DiscreteTrained symbol sh prob+   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 $+          Array.map+             (StorableArray.accumulate (+)+                 (Vector.constant (Shape.ZeroBased $ length set) 0) .+              map (mapFst+                 (mapLookup "estimateDiscrete: unknown emission symbol" emi)))+             grouped+   combine (DiscreteTrained distr0) (DiscreteTrained distr1) =+      DiscreteTrained $ Map.unionWith Vector.add distr0 distr1+   normalize (DiscreteTrained distr) =+      Discrete $ if Map.null distr then distr else normalizeProbVecs distr++transposeVectorList ::+   (Shape.C sh, Eq sh, Class.Real a) =>+   Array sh (Vector Matrix.ZeroInt a) -> [Vector sh a]+transposeVectorList xs =+   case Array.toList xs of+      [] -> []+      x:_ -> Matrix.toRows $ Matrix.fromColumnArray (StorableArray.shape x) xs++normalizeProbVecs ::+   (Shape.C sh, Eq sh, Foldable f, Functor f, Class.Real a) =>+   f (Vector sh a) -> f (Vector sh a)+normalizeProbVecs vs =+   let factors =+         StorableArray.map recip $ List.foldl1' Vector.add $ Fold.toList vs+   in fmap (Vector.mul factors) vs++mapLookup :: (Ord k) => String -> Map.Map k a -> k -> a+mapLookup msg dict x = Map.findWithDefault (error msg) x dict+++newtype Gaussian emiSh stateSh a =+      Gaussian (Array stateSh (Vector emiSh a, UpperTriangular emiSh a, a))+   deriving (Show)++newtype GaussianTrained emiSh stateSh a =+   GaussianTrained+      (Array stateSh+         (Maybe (Vector emiSh a, HermitianMatrix emiSh a, a)))+   deriving (Show)++type instance Probability (Gaussian emiSh stateSh a) = a+type instance Emission (Gaussian emiSh stateSh a) = Vector emiSh a+type instance StateShape (Gaussian emiSh stateSh a) = stateSh+++instance+   (NFData emiSh, NFData stateSh, Shape.C stateSh, NFData a, Storable a) =>+      NFData (Gaussian emiSh stateSh a) where+   rnf (Gaussian params) = rnf params++instance+   (NFData emiSh, NFData stateSh, Shape.C stateSh, NFData a, Storable a) =>+      NFData (GaussianTrained emiSh stateSh a) where+   rnf (GaussianTrained params) = rnf params+++instance+   (FormatArray emiSh, Shape.C stateSh, Class.Real a) =>+      Format (Gaussian emiSh stateSh a) where+   format = runFormatGaussian $ Class.switchReal formatGaussian formatGaussian++newtype FormatGaussian emiSh stateSh a =+   FormatGaussian+      {runFormatGaussian :: String -> Gaussian emiSh stateSh a -> TextBox.Box}++formatGaussian ::+   (FormatArray emiSh, Shape.C stateSh, Class.Real a, Format a) =>+   FormatGaussian emiSh stateSh a+formatGaussian =+   FormatGaussian $ \fmt (Gaussian params) -> format fmt $ Array.toList params+++instance+   (Shape.Indexed stateSh, Eq stateSh, Class.Real a) =>+      Info (Gaussian emiSh stateSh a) where+   statesShape (Gaussian params) = Array.shape params++instance+   (Shape.C emiSh, Eq emiSh, Shape.Indexed stateSh, Eq stateSh, Class.Real a) =>+      Generate (Gaussian emiSh stateSh a) where+   generate (Gaussian allParams) state = do+      let (center, covarianceChol, _c) = allParams ! state+      seed <- MS.state Rnd.random+      return $+         Vector.add center $+         Vector.random Vector.Normal (StorableArray.shape center) seed+            <# covarianceChol++instance+   (Shape.C emiSh, Eq emiSh, Shape.Indexed stateSh, Eq stateSh, Class.Real a) =>+      EmissionProb (Gaussian emiSh stateSh a) where+   emissionProb (Gaussian allParams) x =+      Vector.fromList (Array.shape allParams) $+      map (emissionProbGen x) $ Array.toList allParams+   emissionStateProb (Gaussian allParams) x s =+      emissionProbGen x $ allParams ! s++emissionProbGen ::+   (Shape.C emiSh, Eq emiSh, Class.Real a) =>+   Vector emiSh a -> (Vector emiSh a, UpperTriangular emiSh a, a) -> a+emissionProbGen x (center, covarianceChol, c) =+   let x0 =+         Matrix.solveVector (Triangular.transpose covarianceChol) $+         Vector.sub x center+   in  c * cexp ((-1/2) * Vector.inner x0 x0)+++instance+   (Shape.C emiSh, Eq emiSh, Shape.Indexed stateSh, Eq stateSh, Class.Real a) =>+      Estimate+         (GaussianTrained emiSh stateSh a)+         (Gaussian emiSh stateSh a) where+   type Distribution (GaussianTrained emiSh stateSh a) =+            Gaussian emiSh stateSh a+   type Trained (Gaussian emiSh stateSh a) = GaussianTrained emiSh stateSh a+   accumulateEmissions =+      let params xs =+            (NonEmpty.foldl1Map Vector.add (uncurry $ flip Vector.scale) xs,+             covarianceReal $ fmap swap xs,+             Fold.sum $ fmap snd xs)+      in  GaussianTrained . fmap (fmap params . NonEmpty.fetch)+   combine (GaussianTrained distr0) (GaussianTrained distr1) =+      let comb (center0, covariance0, weight0)+               (center1, covariance1, weight1) =+             (Vector.add center0 center1,+              Vector.add covariance0 covariance1,+              weight0 + weight1)+      in  GaussianTrained $ Array.zipWith (maybePlus comb) distr0 distr1+   {-+     Sum_i (xi-m) * (xi-m)^T+   = Sum_i xi*xi^T + Sum_i m*m^T - Sum_i xi*m^T - Sum_i m*xi^T+   = Sum_i xi*xi^T - Sum_i m*m^T+   = Sum_i xi*xi^T - n * m*m^T+   -}+   normalize (GaussianTrained distr) =+      let params (centerSum, covarianceSum, weight) =+             let c = recip weight+                 center = Vector.scale c centerSum+             in  (center,+                  Vector.sub (Vector.scale c covarianceSum)+                     (Hermitian.outer MatrixShape.RowMajor center))+      in  Gaussian $+          fmap+             (gaussianParameters . params .+              fromMaybe+                (error "Distribution.normalize: undefined array element")) $+          distr++-- ToDo: could be managed by semigroup+maybePlus :: (a -> a -> a) -> Maybe a -> Maybe a -> Maybe a+maybePlus f mx my = liftA2 f mx my <|> mx <|> my+++newtype CovarianceReal f emiSh a =+   CovarianceReal+      {getCovarianceReal :: f (a, Vector emiSh a) -> HermitianMatrix emiSh a}++covarianceReal ::+   (Shape.C emiSh, Eq emiSh, Class.Real a) =>+   NonEmpty.T [] (a, Vector emiSh a) -> HermitianMatrix emiSh a+covarianceReal =+   getCovarianceReal $+   Class.switchReal+      (CovarianceReal $ Hermitian.sumRank1NonEmpty MatrixShape.RowMajor)+      (CovarianceReal $ Hermitian.sumRank1NonEmpty MatrixShape.RowMajor)++gaussian ::+   (Shape.C emiSh, Shape.C stateSh, Class.Real prob) =>+   Array stateSh (Vector emiSh prob, HermitianMatrix emiSh prob) ->+   Gaussian emiSh stateSh prob+gaussian = consGaussian . fmap gaussianParameters++gaussianParameters ::+   (Shape.C emiSh, Class.Real prob) =>+   (Vector emiSh prob, HermitianMatrix emiSh prob) ->+   (Vector emiSh prob, UpperTriangular emiSh prob, prob)+gaussianParameters (center, covariance) =+   gaussianFromCholesky center $ HermitianPD.decompose covariance++consGaussian ::+   (Shape.C stateSh) =>+   Array stateSh (Vector emiSh a, UpperTriangular emiSh a, a) ->+   Gaussian emiSh stateSh a+consGaussian = Gaussian++gaussianFromCholesky ::+   (Shape.C emiSh, Class.Real prob) =>+   Vector emiSh prob -> UpperTriangular emiSh prob ->+   (Vector emiSh prob, UpperTriangular emiSh prob, prob)+gaussianFromCholesky center covarianceChol =+   let covarianceSqrtDet =+         Vector.product $ Triangular.takeDiagonal covarianceChol+   in  (center, covarianceChol,+        recip (sqrt2pi ^ vectorDim center * covarianceSqrtDet))++sqrt2pi :: (Class.Real a) => a+sqrt2pi = runIdentity $ Class.switchReal sqrt2piAux sqrt2piAux++sqrt2piAux :: (Floating a) => Identity a+sqrt2piAux = Identity $ sqrt (2*pi)++cexp :: (Class.Real a) => a -> a+cexp = appEndo $ Class.switchReal (Endo exp) (Endo exp)++++class ToCSV distr where+   toCells :: distr -> [[String]]++class FromCSV distr where+   parseCells :: StateShape distr -> 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+   (Shape.C sh, Class.Real prob, Show prob, Read prob, CSVSymbol symbol) =>+      ToCSV (Discrete symbol sh prob) where+   toCells (Discrete m) =+      map+         (\(symbol, probs) ->+            cellFromSymbol symbol : HMMCSV.cellsFromVector probs) $+      Map.toAscList m++instance+   (Shape.C sh, Class.Real prob, Show prob, Read prob, CSVSymbol symbol) =>+      FromCSV (Discrete symbol sh prob) where+   parseCells n =+      fmap (Discrete . Map.fromList) $+      HMMCSV.manyRowsUntilEnd $ parseSymbolProb n++parseSymbolProb ::+   (Shape.C sh, Class.Real prob, Read prob, CSVSymbol symbol) =>+   sh -> CSV.CSVRow -> HMMCSV.CSVParser (symbol, Vector sh prob)+parseSymbolProb sh 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+                let n = Shape.size sh+                let m = vectorDim v+                HMMCSV.assert (n == m)+                   (printf "number of states (%d) and size of probability vector (%d) mismatch"+                      n m)+                return $ StorableArray.reshape sh v)+++instance+   (Shape.Indexed emiSh, Shape.Indexed stateSh,+    Class.Real a, Eq a, Show a, Read a) =>+      ToCSV (Gaussian emiSh stateSh a) where+   toCells (Gaussian params) =+      List.intercalate [[]] $+      map+         (\(center, covarianceChol, _) ->+            HMMCSV.cellsFromVector center :+            HMMCSV.cellsFromSquare (Triangular.toSquare covarianceChol)) $+      Array.toList params++instance+   (emiSh ~ Matrix.ZeroInt, Shape.Indexed stateSh,+    Class.Real a, Eq a, Show a, Read a) =>+      FromCSV (Gaussian emiSh stateSh a) where+   parseCells sh = do+      let n = Shape.size sh+      gs <- HMMCSV.manySepUntilEnd parseSingleGaussian+      HMMCSV.assert (length gs == n) $+         printf "number of states (%d) and number of Gaussians (%d) mismatch"+            n (length gs)+      let sizes = map (vectorDim . fst3) gs+      HMMCSV.assert (ListHT.allEqual sizes) $+         printf "dimensions of emissions mismatch: %s" (show sizes)+      return $ consGaussian $ Array.fromList sh gs++parseSingleGaussian ::+   (emiSh ~ Matrix.ZeroInt, Class.Real prob, Eq prob, Read prob) =>+   HMMCSV.CSVParser (Vector emiSh prob, UpperTriangular emiSh prob, prob)+parseSingleGaussian = do+   center <- HMMCSV.parseNonEmptyVectorCells+   covarianceCholSquare <-+      HMMCSV.parseSquareMatrixCells $ StorableArray.shape center+   let covarianceChol = Triangular.takeUpper covarianceCholSquare+   HMMCSV.assert+      (isUpperTriang covarianceCholSquare 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 ::+   (Shape.C sh, Class.Real a, Eq a) =>+   SquareMatrix sh a -> UpperTriangular sh a -> Bool+isUpperTriang m mt =+   Vector.toList m == Vector.toList (Triangular.toSquare mt)

@@ -0,0 +1,12 @@+{- |+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" #-}+   (module Math.HiddenMarkovModel.Example.CirclePrivate) where++import Math.HiddenMarkovModel.Example.CirclePrivate

@@ -0,0 +1,123 @@+module Math.HiddenMarkovModel.Example.CirclePrivate where++import qualified Math.HiddenMarkovModel as HMM+import qualified Math.HiddenMarkovModel.Distribution as Distr+import Math.HiddenMarkovModel.Utility+         (normalizeProb, squareFromLists, hermitianFromList)++import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Vector (Vector)++import qualified Data.Array.Comfort.Boxed as Array+import qualified Data.Array.Comfort.Shape as Shape++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 ((!:))+import Data.Maybe (fromMaybe)++++data State = Q1 | Q2 | Q3 | Q4+   deriving (Eq, Ord, Enum, Bounded)++type StateSet = Shape.Enumeration State++stateSet :: StateSet+stateSet = Shape.Enumeration+++data Coordinate = X | Y+   deriving (Eq, Ord, Enum, Bounded)++type CoordinateSet = Shape.Enumeration Coordinate++coordinateSet :: CoordinateSet+coordinateSet = Shape.Enumeration++type HMM = HMM.Gaussian CoordinateSet StateSet Double++hmm :: HMM+hmm =+   HMM.Cons {+      HMM.initial = normalizeProb $ Vector.constant stateSet 1,+      HMM.transition =+         squareFromLists stateSet $+            stateVector 0.9 0.0 0.0 0.1 :+            stateVector 0.1 0.9 0.0 0.0 :+            stateVector 0.0 0.1 0.9 0.0 :+            stateVector 0.0 0.0 0.1 0.9 :+            [],+      HMM.distribution =+         let cov0 = hermitianFromList coordinateSet [0.10, -0.09, 0.10]+             cov1 = hermitianFromList coordinateSet [0.10,  0.09, 0.10]+         in  Distr.gaussian $ Array.fromList stateSet $+                (Vector.fromList coordinateSet [ 0.5,  0.5], cov0) :+                (Vector.fromList coordinateSet [-0.5,  0.5], cov1) :+                (Vector.fromList coordinateSet [-0.5, -0.5], cov0) :+                (Vector.fromList coordinateSet [ 0.5, -0.5], cov1) :+                []+   }++stateVector :: Double -> Double -> Double -> Double -> Vector StateSet Double+stateVector x0 x1 x2 x3 = Vector.fromList stateSet [x0,x1,x2,x3]++circleLabeled :: NonEmpty.T [] (State, Vector CoordinateSet Double)+circleLabeled =+   NonEmpty.mapTail (take 200) $+   fmap+      (\x ->+         (toEnum $ mod (floor (x*2/pi)) 4,+          Vector.fromList coordinateSet [cos x, sin x])) $+   NonEmptyC.iterate (0.1+) 0++circle :: NonEmpty.T [] (Vector CoordinateSet Double)+circle = fmap snd circleLabeled++revealed :: NonEmpty.T [] 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 CoordinateSet () Double+reconstructDistribution =+   let gen = Distr.generate (HMM.distribution hmm) Q1+   in  HMM.finishTraining $ HMM.trainSupervised () $ fmap ((,) ()) $+       flip MS.evalState (Rnd.mkStdGen 23) $+       liftM2 (!:) gen $ replicateM 1000 gen++{- |+Generate labeled emission sequences+and use them for supervised training.+-}+reconstructModel :: HMM+reconstructModel =+   HMM.trainMany (HMM.trainSupervised stateSet) $+   fmap+      (\seed ->+         fromMaybe (error "empty generated sequence") $ NonEmpty.fetch $+         take 1000 $ HMM.generateLabeled hmm $ Rnd.mkStdGen seed)+      (23 !: take 42 [24..])+++hmmTrainedSupervised :: HMM+hmmTrainedSupervised =+   HMM.finishTraining $ HMM.trainSupervised stateSet circleLabeled++hmmTrainedUnsupervised :: HMM+hmmTrainedUnsupervised =+   HMM.finishTraining $ HMM.trainUnsupervised hmm circle++hmmIterativelyTrained :: HMM+hmmIterativelyTrained =+   nest 100+      (HMM.finishTraining . flip HMM.trainUnsupervised circle)+      hmm

@@ -0,0 +1,90 @@+{- |+Example of an HMM with continuous emissions.+We train a model to accept sine waves of a certain frequency.++There are four hidden states: 'Rising', 'High', 'Falling', '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 Math.HiddenMarkovModel.Utility+         (normalizeProb, squareFromLists, hermitianFromList, singleton)++import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Vector (Vector)++import qualified Data.Array.Comfort.Boxed as Array+import qualified Data.Array.Comfort.Shape as Shape++import qualified Data.NonEmpty.Class as NonEmptyC+import qualified Data.NonEmpty as NonEmpty+import Data.Function.HT (nest)+import Data.Tuple.HT (mapSnd)++++data State = Rising | High | Falling | Low+   deriving (Eq, Ord, Enum, Bounded)++type StateSet = Shape.Enumeration State++stateSet :: StateSet+stateSet = Shape.Enumeration+++type HMM = HMM.Gaussian () StateSet Double++hmm :: HMM+hmm =+   HMM.Cons {+      HMM.initial = normalizeProb $ Vector.constant stateSet 1,+      HMM.transition =+         squareFromLists stateSet $+            stateVector 0.9 0.0 0.0 0.1 :+            stateVector 0.1 0.9 0.0 0.0 :+            stateVector 0.0 0.1 0.9 0.0 :+            stateVector 0.0 0.0 0.1 0.9 :+            [],+      HMM.distribution =+         Distr.gaussian $ Array.fromList stateSet $+            (singleton   0 , hermitianFromList () [1]) :+            (singleton   1 , hermitianFromList () [1]) :+            (singleton   0 , hermitianFromList () [1]) :+            (singleton (-1), hermitianFromList () [1]) :+            []+   }++stateVector :: Double -> Double -> Double -> Double -> Vector StateSet Double+stateVector x0 x1 x2 x3 = Vector.fromList stateSet [x0,x1,x2,x3]++sineWaveLabeled :: NonEmpty.T [] (State, Double)+sineWaveLabeled =+   NonEmpty.mapTail (take 200) $+   fmap (\x -> (toEnum $ 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 [] State+revealed = HMM.reveal hmmTrainedSupervised $ fmap singleton sineWave++hmmTrainedSupervised :: HMM+hmmTrainedSupervised =+   HMM.finishTraining $ HMM.trainSupervised stateSet $+   fmap (mapSnd singleton) sineWaveLabeled++hmmTrainedUnsupervised :: HMM+hmmTrainedUnsupervised =+   HMM.finishTraining $ HMM.trainUnsupervised hmm $ fmap singleton sineWave++hmmIterativelyTrained :: HMM+hmmIterativelyTrained =+   nest 100+      (\model ->+         HMM.finishTraining $ HMM.trainUnsupervised model $+         fmap singleton sineWave)+      hmm

@@ -0,0 +1,50 @@+{- |+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:+'StateRed' emits red, 'StateYellowRG' emits yellow between red and green,+'StateGreen' emits green, 'StateYellowGR' 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" #-}+   (+   HMM,+   Color(..),+   hmm,+   hmmDisturbed,+   red, yellowRG, green, yellowGR,+   labeledSequences,+   hmmTrainedSupervised,+   stateSequences,+   hmmTrainedUnsupervised,+   hmmIterativelyTrained,+   ) where++import Math.HiddenMarkovModel.Example.TrafficLightPrivate

@@ -0,0 +1,164 @@+{-# LANGUAGE TypeFamilies #-}+module Math.HiddenMarkovModel.Example.TrafficLightPrivate where++import qualified Math.HiddenMarkovModel as HMM+import qualified Math.HiddenMarkovModel.Distribution as Distr+import Math.HiddenMarkovModel.Utility (normalizeProb, squareFromLists)++import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Format (Format(format))++import qualified Data.Array.Comfort.Shape as Shape++import qualified Text.PrettyPrint.Boxes as TextBox+import Text.Read.HT (maybeRead)++import Control.DeepSeq (NFData(rnf))+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)++instance NFData Color where+   rnf Red = ()+   rnf _ = ()++instance Format Color where+   format _fmt = TextBox.text . show++{- |+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+++-- data State = StateRed | StateYellowDown | StateGreen | StateYellowUp+data State = StateRed | StateYellowRG | StateGreen | StateYellowGR+   deriving (Eq, Ord, Enum, Bounded)++type StateSet = Shape.Enumeration State++stateSet :: StateSet+stateSet = Shape.Enumeration+++type HMM = HMM.Discrete Color StateSet Double++hmm :: HMM+hmm =+   HMM.Cons {+      HMM.initial = normalizeProb $ stateVector 2 1 2 1,+      HMM.transition =+         squareFromLists stateSet $+            stateVector 0.8 0.0 0.0 0.2 :+            stateVector 0.2 0.8 0.0 0.0 :+            stateVector 0.0 0.2 0.8 0.0 :+            stateVector 0.0 0.0 0.2 0.8 :+            [],+      HMM.distribution =+         Distr.Discrete $ Map.fromList $+            (Red,    stateVector 1 0 0 0) :+            (Yellow, stateVector 0 1 0 1):+            (Green,  stateVector 0 0 1 0) :+            []+   }++hmmDisturbed :: HMM+hmmDisturbed =+   HMM.Cons {+      HMM.initial = normalizeProb $ stateVector 1 1 1 1,+      HMM.transition =+         squareFromLists stateSet $+            stateVector 0.3 0.2 0.2 0.3 :+            stateVector 0.3 0.3 0.2 0.2 :+            stateVector 0.2 0.3 0.3 0.2 :+            stateVector 0.2 0.2 0.3 0.3 :+            [],+      HMM.distribution =+         Distr.Discrete $ Map.fromList $+            (Red,    stateVector 0.6 0.2 0.2 0.2) :+            (Yellow, stateVector 0.2 0.6 0.2 0.6) :+            (Green,  stateVector 0.2 0.2 0.6 0.2) :+            []+   }++stateVector :: Double -> Double -> Double -> Double -> Vector StateSet Double+stateVector x0 x1 x2 x3 = Vector.fromList stateSet [x0,x1,x2,x3]+++red, yellowRG, green, yellowGR :: (State, Color)+red      = (StateRed, Red)+yellowRG = (StateYellowRG, Yellow)+green    = (StateGreen, Green)+yellowGR = (StateYellowGR, Yellow)++labeledSequences :: NonEmpty.T [] (NonEmpty.T [] (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+hmmTrainedSupervised =+   HMM.trainMany (HMM.trainSupervised stateSet) 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+hmmTrainedUnsupervised =+   HMM.trainMany (HMM.trainUnsupervised hmm) stateSequences++{- |+Repeat unsupervised training until convergence.+-}+hmmIterativelyTrained :: HMM+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 -> NonEmpty.T [] (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)

@@ -0,0 +1,117 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+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.Utility (attachOnes, vectorDim)++import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable as StorableArray+import qualified Data.Array.Comfort.Boxed as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Boxed (Array)++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 Control.DeepSeq (NFData, rnf)+import Foreign.Storable (Storable)++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 sh ix prob =+   Cons {+      model :: HMM.T distr sh prob,+      nameFromStateMap :: Array sh String,+      stateFromNameMap :: Map String ix+   }+   deriving (Show)++type Discrete symbol stateSh prob =+      T (Distr.Discrete symbol stateSh prob) stateSh (Shape.Index stateSh) prob+type Gaussian emiSh stateSh a =+      T (Distr.Gaussian emiSh stateSh a) stateSh (Shape.Index stateSh) a+++instance+   (NFData distr, NFData sh, NFData ix, NFData prob,+    Shape.C sh, Storable prob) =>+      NFData (T distr sh ix prob) where+   rnf hmm = rnf (model hmm, nameFromStateMap hmm, stateFromNameMap hmm)+++fromModelAndNames ::+   (Shape.Indexed sh, Shape.Index sh ~ state) =>+   HMM.T distr sh prob -> [String] -> T distr sh state prob+fromModelAndNames md names =+   let m = Array.fromList (StorableArray.shape $ HMM.initial md) names+   in  Cons {+          model = md,+          nameFromStateMap = m,+          stateFromNameMap = inverseMap m+       }++inverseMap ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) => Array sh String -> Map String ix+inverseMap =+   Map.fromListWith (error "duplicate label") .+   map swap . Array.toAssociations+++toCSV ::+   (Distr.ToCSV distr, Shape.Indexed sh, Class.Real prob, Show prob) =>+   T distr sh ix prob -> String+toCSV hmm =+   CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $+      Array.toList (nameFromStateMap hmm) : HMM.toCells (model hmm)++fromCSV ::+   (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh,+    Shape.Indexed stateSh, Shape.Index stateSh ~ state,+    Class.Real prob, Read prob) =>+   (Int -> stateSh) ->+   String -> ME.Exceptional String (T distr stateSh state prob)+fromCSV makeShape =+   MS.evalStateT (parseCSV makeShape) . map HMMCSV.fixShortRow . CSV.parseCSV++parseCSV ::+   (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh,+    Shape.Indexed stateSh, Shape.Index stateSh ~ state,+    Class.Real prob, Read prob) =>+   (Int -> stateSh) -> HMMCSV.CSVParser (T distr stateSh state prob)+parseCSV makeShape = 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 makeShape+   let n = length names+       m = vectorDim (HMM.initial md)+    in HMMCSV.assert (n == m) $+          printf "got %d state names for %d states" n m+   return $ fromModelAndNames md names

@@ -0,0 +1,175 @@+{-# LANGUAGE TypeFamilies #-}+{- |+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,+           biscaleTransition, matrixMaxMul, sumTransitions)+import Math.HiddenMarkovModel.Utility+         (SquareMatrix, normalizeFactor, normalizeProb)++import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix ((<#), (#>))+import Numeric.LAPACK.Vector (Vector)++import qualified Numeric.Netlib.Class as Class++import qualified Control.Functor.HT as Functor++import qualified Data.Array.Comfort.Storable as StorableArray+import qualified Data.Array.Comfort.Boxed as Array+import qualified Data.Array.Comfort.Shape as Shape++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, Distr.StateShape distr ~ sh, Eq sh, Floating prob,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr sh prob -> NonEmpty.T f emission -> prob+logLikelihood hmm = Fold.sum . fmap (log . fst) . alpha hmm++alpha ::+   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr sh prob ->+   NonEmpty.T f emission -> NonEmpty.T f (prob, Vector sh prob)+alpha hmm (NonEmpty.Cons x xs) =+   let normMulEmiss y = normalizeFactor . Vector.mul (emission hmm y)+   in  NonEmpty.scanl+          (\(_,alphai) xi -> normMulEmiss xi (transition hmm #> alphai))+          (normMulEmiss x (initial hmm))+          xs++beta ::+   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f, NonEmptyC.Reverse f) =>+   T distr sh prob ->+   f (prob, emission) -> NonEmpty.T f (Vector sh prob)+beta hmm =+   nonEmptyScanr+      (\(ci,xi) betai ->+         Vector.scale (recip ci) $+         Vector.mul (emission hmm xi) betai <# transition hmm)+      (Vector.constant (StorableArray.shape $ initial hmm) 1)++alphaBeta ::+   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f, NonEmptyC.Zip f, NonEmptyC.Reverse f) =>+   T distr sh prob ->+   NonEmpty.T f emission ->+   (NonEmpty.T f (prob, Vector sh prob), NonEmpty.T f (Vector sh 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.StateShape distr ~ sh, Eq sh,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f, NonEmptyC.Zip f) =>+   T distr sh prob ->+   NonEmpty.T f emission ->+   NonEmpty.T f (prob, Vector sh prob) ->+   NonEmpty.T f (Vector sh prob) ->+   f (SquareMatrix sh prob)+xiFromAlphaBeta hmm xs calphas betas =+   let (cs,alphas) = Functor.unzip calphas+   in  NonEmptyC.zipWith4+          (\x alpha0 c1 beta1 ->+             Vector.scale (recip c1) $ biscaleTransition hmm x alpha0 beta1)+          (NonEmpty.tail xs)+          (NonEmpty.init alphas)+          (NonEmpty.tail cs)+          (NonEmpty.tail betas)++zetaFromAlphaBeta ::+   (Shape.C sh, Eq sh, Class.Real prob, NonEmptyC.Zip f) =>+   NonEmpty.T f (prob, Vector sh prob) ->+   NonEmpty.T f (Vector sh prob) ->+   NonEmpty.T f (Vector sh prob)+zetaFromAlphaBeta calphas betas =+   NonEmptyC.zipWith (Vector.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.StateShape distr ~ sh,+    Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ state,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f, NonEmptyC.Reverse f) =>+   T distr sh prob -> NonEmpty.T f emission -> NonEmpty.T f state+reveal hmm (NonEmpty.Cons x xs) =+   fmap (Shape.revealIndex (StorableArray.shape $ initial hmm)) $+   uncurry (NonEmpty.scanr (StorableArray.!)) $+   mapFst+      (fst . Vector.argAbsMaximum .+       StorableArray.mapShape Shape.Deferred) $+   mapAccumL+      (\alphai xi ->+         swap $ mapSnd (Vector.mul (emission hmm xi)) $+         matrixMaxMul (transition hmm) $ normalizeProb alphai)+      (Vector.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, Distr.StateShape distr ~ sh, Eq sh,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   T distr sh prob -> NonEmpty.T [] emission -> Trained tdistr sh prob+trainUnsupervised hmm xs =+   let (alphas, betas) = alphaBeta hmm xs+       zetas = zetaFromAlphaBeta alphas betas+       zeta0 = NonEmpty.head zetas++   in  Trained {+          trainedInitial = zeta0,+          trainedTransition =+             sumTransitions hmm $ xiFromAlphaBeta hmm xs alphas betas,+          trainedDistribution =+             Distr.accumulateEmissions $+             Array.fromList (StorableArray.shape zeta0) $+             map (zip (NonEmpty.flatten xs)) $+             List.transpose $ map Vector.toList $ NonEmpty.flatten zetas+       }

@@ -0,0 +1,121 @@+{-# LANGUAGE TypeFamilies #-}+{- |+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:++> example :: HMM.DiscreteTrained Char (ShapeStatic.ZeroBased TypeNum.U2) Double+> example =+>    let a = atom FL.i0+>        b = atom FL.i1+>        distr =+>           Distr.DiscreteTrained $ Map.fromList $+>           ('a', ShapeStatic.vector $ 1!:2!:FL.end) :+>           ('b', ShapeStatic.vector $ 4!:3!:FL.end) :+>           ('c', ShapeStatic.vector $ 0!:1!:FL.end) :+>           []+>    in finish (ShapeStatic.ZeroBased Proxy) 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.Utility (SquareMatrix, squareConstant)++import qualified Numeric.LAPACK.Vector as Vector+import qualified Numeric.LAPACK.ShapeStatic as ShapeStatic++import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable as StorableArray+import qualified Data.Array.Comfort.Shape as Shape++import qualified Data.FixedLength as FL+import Data.FixedLength ((!:))++import qualified Type.Data.Num.Unary.Literal as TypeNum+import Type.Base.Proxy (Proxy(Proxy))++import qualified Data.Map as Map+import Data.Semigroup (Semigroup, (<>), stimes)++import Prelude hiding (replicate)+++newtype T sh prob =+   Cons (sh -> (Shape.Index sh, SquareMatrix sh prob, Shape.Index sh))++atom ::+   (Shape.Indexed sh, Shape.Index sh ~ state, Class.Real prob) =>+   state -> T sh prob+atom s = Cons $ \sh -> (s, squareConstant sh 0, s)+++instance+   (Shape.Indexed sh, Eq sh, Class.Real prob) =>+      Semigroup (T sh prob) where+   (<>) = append+   stimes k = replicate $ fromIntegral k+++infixl 5 `append`++append ::+   (Shape.Indexed sh, Eq sh, Class.Real prob) =>+   T sh prob -> T sh prob -> T sh 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 $ Vector.add ma mb, sbo)++replicate ::+   (Shape.Indexed sh, Class.Real prob) => Int -> T sh prob -> T sh prob+replicate ki (Cons f) =+   Cons $ \sh ->+      case f sh of+         (si, m, so) ->+            let k = fromIntegral ki+            in  (si, increment (si,so) (k-1) $ Vector.scale k m, so)++increment ::+   (Shape.Indexed sh, Shape.Index sh ~ state, Class.Real a) =>+   (state, state) -> a -> SquareMatrix sh a -> SquareMatrix sh a+increment (i,j) x m  =  StorableArray.accumulate (+) m [((i,j), x)]+++finish ::+   (Shape.Indexed sh, Class.Real prob) =>+   sh -> tdistr -> T sh prob -> Trained tdistr sh prob+finish sh tdistr (Cons f) =+   case f sh of+      (si, m, _so) ->+         Trained {+            trainedInitial = StorableArray.fromAssociations sh 0 [(si,1)],+            trainedTransition = m,+            trainedDistribution = tdistr+         }+++_example :: HMM.DiscreteTrained Char (ShapeStatic.ZeroBased TypeNum.U2) Double+_example =+   let a = atom FL.i0+       b = atom FL.i1+       distr =+          Distr.DiscreteTrained $ Map.fromList $+          ('a', ShapeStatic.vector $ 1!:2!:FL.end) :+          ('b', ShapeStatic.vector $ 4!:3!:FL.end) :+          ('c', ShapeStatic.vector $ 0!:1!:FL.end) :+          []+   in finish (ShapeStatic.ZeroBased Proxy) distr $+      replicate 5 $ replicate 10 a <> replicate 20 b

@@ -0,0 +1,335 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UndecidableInstances #-}+module Math.HiddenMarkovModel.Private where++import qualified Math.HiddenMarkovModel.Distribution as Distr+import qualified Math.HiddenMarkovModel.CSV as HMMCSV+import Math.HiddenMarkovModel.Utility (SquareMatrix, diagonal)++import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import qualified Numeric.LAPACK.Format as Format+import Numeric.LAPACK.Matrix ((<#), (<#>), (#>))+import Numeric.LAPACK.Vector (Vector)++import qualified Numeric.Netlib.Class as Class++import Control.DeepSeq (NFData, rnf)+import Control.Applicative ((<$>))++import Foreign.Storable (Storable)++import qualified Data.Array.Comfort.Storable as StorableArray+import qualified Data.Array.Comfort.Boxed as Array+import qualified Data.Array.Comfort.Shape as Shape++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 sh prob =+   Cons {+      initial :: Vector sh prob,+      transition :: SquareMatrix sh prob,+      distribution :: distr+   }+   deriving (Show)++instance+   (NFData distr, NFData sh, NFData prob, Storable prob) =>+      NFData (T distr sh prob) where+   rnf (Cons initial_ transition_ distribution_) =+      rnf (initial_, transition_, distribution_)++instance+   (Class.Real prob, Format.FormatArray sh, Format.Format distr) =>+      Format.Format (T distr sh prob) where+   format fmt (Cons initial_ transition_ distribution_) =+      Format.format fmt (initial_, transition_, distribution_)+++emission ::+   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   T distr sh prob -> emission -> Vector sh prob+emission  =  Distr.emissionProb . distribution+++forward ::+   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr sh prob -> NonEmpty.T f emission -> prob+forward hmm = Vector.sum . NonEmpty.last . alpha hmm++alpha ::+   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr sh prob ->+   NonEmpty.T f emission -> NonEmpty.T f (Vector sh prob)+alpha hmm (NonEmpty.Cons x xs) =+   NonEmpty.scanl+      (\alphai xi -> Vector.mul (emission hmm xi) (transition hmm #> alphai))+      (Vector.mul (emission hmm x) (initial hmm))+      xs+++backward ::+   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr sh prob -> NonEmpty.T f emission -> prob+backward hmm (NonEmpty.Cons x xs) =+   Vector.sum $+   Vector.mul (initial hmm) $+   Vector.mul (emission hmm x) $+   NonEmpty.head $ beta hmm xs++beta ::+   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr sh prob ->+   f emission -> NonEmpty.T f (Vector sh prob)+beta hmm =+   NonEmpty.scanr+      (\xi betai -> Vector.mul (emission hmm xi) betai <# transition hmm)+      (Vector.constant (StorableArray.shape $ initial hmm) 1)+++alphaBeta ::+   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr sh prob ->+   NonEmpty.T f emission ->+   (prob, NonEmpty.T f (Vector sh prob), NonEmpty.T f (Vector sh prob))+alphaBeta hmm xs =+   let alphas = alpha hmm xs+       betas = beta hmm $ NonEmpty.tail xs+       recipLikelihood = recip $ Vector.sum $ NonEmpty.last alphas+   in  (recipLikelihood, alphas, betas)++++biscaleTransition ::+   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr,+    Distr.EmissionProb distr, Distr.Probability distr ~ prob) =>+   T distr sh prob -> Distr.Emission distr ->+   Vector sh prob -> Vector sh prob -> SquareMatrix sh prob+biscaleTransition hmm x alpha0 beta1 =+   diagonal (Vector.mul (emission hmm x) beta1)+   <#>+   transition hmm+   <#>+   diagonal alpha0++xiFromAlphaBeta ::+   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   T distr sh prob -> prob ->+   NonEmpty.T [] emission ->+   NonEmpty.T [] (Vector sh prob) ->+   NonEmpty.T [] (Vector sh prob) ->+   [SquareMatrix sh prob]+xiFromAlphaBeta hmm recipLikelihood xs alphas betas =+   zipWith3+      (\x alpha0 beta1 ->+         Vector.scale recipLikelihood $+         biscaleTransition hmm x alpha0 beta1)+      (NonEmpty.tail xs)+      (NonEmpty.init alphas)+      (NonEmpty.tail betas)++zetaFromXi ::+   (Shape.C sh, Eq sh, Class.Real prob) =>+   [SquareMatrix sh prob] -> [Vector sh prob]+zetaFromXi = map Matrix.columnSums++zetaFromAlphaBeta ::+   (Shape.C sh, Eq sh, Class.Real prob) =>+   prob ->+   NonEmpty.T [] (Vector sh prob) ->+   NonEmpty.T [] (Vector sh prob) ->+   NonEmpty.T [] (Vector sh prob)+zetaFromAlphaBeta recipLikelihood alphas betas =+   fmap (Vector.scale recipLikelihood) $+   NonEmptyC.zipWith Vector.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 ::+   (Shape.InvIndexed sh, Shape.Index sh ~ state,+    Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+    Traversable f) =>+   T distr sh prob -> NonEmpty.T f emission -> NonEmpty.T f state+reveal hmm (NonEmpty.Cons x xs) =+   fmap (Shape.revealIndex (StorableArray.shape $ initial hmm)) $+   uncurry (NonEmpty.scanr (StorableArray.!)) $+   mapFst+      (fst . Vector.argAbsMaximum .+       StorableArray.mapShape Shape.Deferred) $+   mapAccumL+      (\alphai xi ->+         swap $ mapSnd (Vector.mul (emission hmm xi)) $+         matrixMaxMul (transition hmm) alphai)+      (Vector.mul (emission hmm x) (initial hmm)) xs++matrixMaxMul ::+   (Shape.Indexed sh, Eq sh, Shape.Index sh ~ ix, Class.Real a) =>+   SquareMatrix sh a -> Vector sh a ->+   (Vector (Shape.Deferred sh) (Shape.DeferredIndex ix), Vector sh a)+matrixMaxMul m v =+   let sh = StorableArray.shape v+   in mapPair (Vector.fromList (Shape.Deferred sh), Vector.fromList sh) $+      unzip $+      map (Vector.argAbsMaximum .+           StorableArray.mapShape Shape.Deferred .+           Vector.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 sh prob =+   Trained {+      trainedInitial :: Vector sh prob,+      trainedTransition :: SquareMatrix sh prob,+      trainedDistribution :: distr+   }+   deriving (Show)++instance+   (NFData distr, NFData sh, NFData prob, Storable prob) =>+      NFData (Trained distr sh prob) where+   rnf hmm =+      rnf (trainedInitial hmm, trainedTransition hmm, trainedDistribution hmm)+++sumTransitions ::+   (Shape.C sh, Eq sh, Class.Real e) =>+   T distr sh e -> [SquareMatrix sh e] -> SquareMatrix sh e+sumTransitions hmm =+   List.foldl' Vector.add+      (Vector.constant (StorableArray.shape $ transition hmm) 0)++{- |+Baum-Welch algorithm+-}+trainUnsupervised ::+   (Distr.Estimate tdistr distr, Distr.StateShape distr ~ sh, Eq sh,+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>+   T distr sh prob -> NonEmpty.T [] emission -> Trained tdistr sh prob+trainUnsupervised hmm xs =+   let (recipLikelihood, alphas, betas) = alphaBeta hmm xs+       zetas = zetaFromAlphaBeta recipLikelihood alphas betas+       zeta0 = NonEmpty.head zetas++   in  Trained {+          trainedInitial = zeta0,+          trainedTransition =+             sumTransitions hmm $+             xiFromAlphaBeta hmm recipLikelihood xs alphas betas,+          trainedDistribution =+             Distr.accumulateEmissions $+             Array.fromList (StorableArray.shape zeta0) $+             map (zip (NonEmpty.flatten xs)) $+             List.transpose $ map Vector.toList $ NonEmpty.flatten zetas+       }+++mergeTrained ::+   (Shape.C sh, Eq sh,+    Distr.Estimate tdistr distr, Distr.Probability distr ~ prob) =>+   Trained tdistr sh prob -> Trained tdistr sh prob -> Trained tdistr sh prob+mergeTrained hmm0 hmm1 =+   Trained {+      trainedInitial = Vector.add (trainedInitial hmm0) (trainedInitial hmm1),+      trainedTransition =+         Vector.add (trainedTransition hmm0) (trainedTransition hmm1),+      trainedDistribution =+         Distr.combine+            (trainedDistribution hmm0) (trainedDistribution hmm1)+   }++instance+   (Shape.C sh, Eq sh,+    Distr.Estimate tdistr distr, Distr.Probability distr ~ prob) =>+      Sg.Semigroup (Trained tdistr sh prob) where+   (<>) = mergeTrained+++toCells ::+   (Distr.ToCSV distr, Shape.Indexed sh, Class.Real prob, Show prob) =>+   T distr sh prob -> [[String]]+toCells hmm =+   (HMMCSV.cellsFromVector $ initial hmm) :+   (HMMCSV.cellsFromSquare $ transition hmm) +++   [] :+   (Distr.toCells $ distribution hmm)++parseCSV ::+   (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh, Shape.C stateSh,+    Class.Real prob, Read prob) =>+   (Int -> stateSh) -> HMMCSV.CSVParser (T distr stateSh prob)+parseCSV makeShape = do+   v <-+      StorableArray.mapShape (makeShape . Shape.zeroBasedSize) <$>+      HMMCSV.parseNonEmptyVectorCells+   let sh = StorableArray.shape v+   m <- HMMCSV.parseSquareMatrixCells sh+   HMMCSV.skipEmptyRow+   distr <- Distr.parseCells sh+   return $ Cons {+      initial = v,+      transition = m,+      distribution = distr+   }

@@ -0,0 +1,258 @@+{- |+Do not import this module, it is only intended for testing!+-}+module Math.HiddenMarkovModel.Test (tests) where++import qualified Math.HiddenMarkovModel.Example.TrafficLightPrivate+                                                            as TrafficLight+import qualified Math.HiddenMarkovModel.Example.CirclePrivate as Circle++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 Math.HiddenMarkovModel.Utility (SquareMatrix, squareFromLists)++import qualified Numeric.LAPACK.Vector as Vector+import qualified Numeric.LAPACK.ShapeStatic as ShapeStatic+import Numeric.LAPACK.Vector (Vector)++import qualified Data.Array.Comfort.Shape as Shape++import qualified Data.FixedLength as FL+import Data.FixedLength ((!:))++import qualified Type.Data.Num.Unary.Literal as TypeNum+import Type.Base.Proxy (Proxy(Proxy))++import qualified Test.QuickCheck as QC+import qualified System.Random as Rnd++import Control.DeepSeq (deepseq)++import qualified Data.NonEmpty.Class as NonEmptyC+import qualified Data.NonEmpty as NonEmpty+import qualified Data.Traversable as Trav+import qualified Data.Foldable as Fold+import qualified Data.Map as Map+import Data.Tuple.HT (mapSnd)++import Text.Printf (printf)+++type StateSet = ShapeStatic.ZeroBased TypeNum.U4++hmm :: HMM.Discrete Char StateSet Double+hmm =+   HMM.Cons {+      HMM.initial = stateVector 0.1 0.2 0.3 0.4,+      HMM.transition =+         squareFromLists stateSet $+            stateVector 0.7 0.1 0.0 0.2 :+            stateVector 0.1 0.6 0.1 0.0 :+            stateVector 0.1 0.2 0.7 0.0 :+            stateVector 0.1 0.1 0.2 0.8 :+            [],+      HMM.distribution =+         Distr.Discrete $ Map.fromList $+            ('a', stateVector 1 0 0 0) :+            ('b', stateVector 0 1 0 1) :+            ('c', stateVector 0 0 1 0) :+            []+   }++stateSet :: StateSet+stateSet = ShapeStatic.ZeroBased Proxy++stateVector :: Double -> Double -> Double -> Double -> Vector StateSet Double+stateVector x0 x1 x2 x3 = ShapeStatic.vector $ x0!:x1!:x2!:x3!:FL.end+-- stateVector = FL.curry ShapeStatic.vector+++sequ :: NonEmpty.T [] Char+sequ = NonEmpty.cons 'a' $ take 20 (HMM.generate hmm (Rnd.mkStdGen 42))++possibleStates :: Char -> [FL.Index TypeNum.U4]+possibleStates c =+   map fst $ filter snd $+   zip (Shape.indices stateSet) $+   map+      (\p ->+         case p of+            0 -> False+            1 -> True+            _ -> error "invalid emission probability (must be 0 or 1)") $+   Vector.toList $+   Map.findWithDefault (error "invalid character") c $+   case HMM.distribution hmm of Distr.Discrete m -> m++{- |+Should all be equal.+-}+sequLikelihood :: ((Double, Double), Double, Double, NonEmpty.T [] Double)+sequLikelihood =+   ((Priv.forward hmm sequ, Priv.backward hmm sequ),+    exp $ Normalized.logLikelihood hmm sequ,+    sum $+       map (NonEmpty.product . HMM.probabilitySequence hmm) $+       Trav.mapM (\c -> map (flip (,) c) $ possibleStates c) sequ,+    NonEmptyC.zipWith Vector.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 Vector.dot (fmap snd calphas) betas+++{- |+Lists should be equal, but the first list contains one less element.+-}+zetas ::+   ([Vector StateSet Double],+    NonEmpty.T [] (Vector StateSet Double),+    NonEmpty.T [] (Vector StateSet Double))+zetas =+   let (recipLikelihood, alphas, betas) = Priv.alphaBeta hmm sequ+   in  (Priv.zetaFromXi $+           Priv.xiFromAlphaBeta hmm recipLikelihood sequ alphas betas,+        Priv.zetaFromAlphaBeta recipLikelihood alphas betas,+        uncurry Normalized.zetaFromAlphaBeta $+        Normalized.alphaBeta hmm sequ)+++distance ::+   (Shape.C sh, Eq sh) =>+   Vector.Vector sh Double -> Vector.Vector sh Double -> Double+distance x y = Vector.normInf (Vector.sub x y)+++{- |+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 $ zipWith distance z0 $ NonEmpty.init z1,+          NonEmpty.maximum $ NonEmptyC.zipWith distance z1 z2)++{- |+Lists should be equal+-}+xis :: ([SquareMatrix StateSet Double], [SquareMatrix StateSet 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 $ zipWith distance x0 x1)+++reveal :: Bool+reveal =+   Normalized.reveal hmm sequ == Priv.reveal hmm sequ+++trainUnsupervised ::+   (HMM.DiscreteTrained Char StateSet Double,+    HMM.DiscreteTrained Char StateSet Double)+trainUnsupervised =+   (Priv.trainUnsupervised hmm sequ,+    Normalized.trainUnsupervised hmm sequ)++trainUnsupervisedDiff :: (Double, Double, (Bool, Double))+trainUnsupervisedDiff =+   case trainUnsupervised of+      (hmm0,hmm1) ->+         (distance (Priv.trainedTransition hmm0) (Priv.trainedTransition hmm1),+          distance+             (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 $ Map.intersectionWith distance m0 m1))+++nonEmptyScanr :: Int -> [Int] -> Bool+nonEmptyScanr x xs =+   Normalized.nonEmptyScanr (-) x xs == NonEmpty.scanr (-) x xs+++circleTraining :: (Int, Circle.HMM) -> Bool+circleTraining (maxDiff,hmm_) =+   maxDiff >=+   (length $ filter id $ NonEmpty.flatten $+    NonEmpty.zipWith (/=)+      (HMM.reveal hmm_ Circle.circle) (fmap fst Circle.circleLabeled))+++allPair :: (a -> Bool, b -> Bool) -> (a,b) -> Bool+allPair (f,g) (a,b) = f a && g b++allTriple :: (a -> Bool, b -> Bool, c -> Bool) -> (a,b,c) -> Bool+allTriple (f,g,h) (a,b,c) = f a && g b && h c++almostZero :: Double -> Bool+almostZero x  =  x < 1e-10++almostOne :: Double -> Bool+almostOne x  =  almostZero $ abs (x-1)++almostEqual :: Double -> Double -> Bool+almostEqual x y  =  almostZero $ abs (x-y)++tests :: [(String, QC.Property)]+tests =+   ("sequLikelihood",+      QC.property $+      case sequLikelihood of+         (forwardBackward, expLog, sumProb, alphaBetas) ->+            allPair (almostEqual sumProb, almostEqual sumProb) forwardBackward+            &&+            almostEqual sumProb expLog+            &&+            length (NonEmpty.tail sequ) == length (NonEmpty.tail alphaBetas)+            &&+            Fold.all (almostEqual sumProb) alphaBetas) :+   ("sequLikelihoodNormalized",+      QC.property $+      length (NonEmpty.tail sequ) ==+         length (NonEmpty.tail sequLikelihoodNormalized)+      &&+      Fold.all almostOne sequLikelihoodNormalized) :+   ("zetasDiff",+      QC.property $ allTriple (id, almostZero, almostZero) zetasDiff) :+   ("xisDiff", QC.property $ allPair (id, almostZero) xisDiff) :+   ("reveal", QC.property reveal) :+   ("trainUnsupervisedDiff",+      QC.property $+      allTriple (almostZero, almostZero, allPair (id, almostZero)) $+      trainUnsupervisedDiff) :+   ("nonEmptyScanr", QC.property nonEmptyScanr) :+   (zip+      (map (printf "TrafficLight.verifyRevelation.%d") [(0::Int) ..])+      (map QC.property TrafficLight.verifyRevelations)) +++   ("TrafficLight.hmmIterativelyTrained.defined",+      QC.property $ deepseq TrafficLight.hmmIterativelyTrained True) :+   (map (mapSnd (QC.property . circleTraining)) $+      ("Circle.hmm", (0, Circle.hmm)) :+      ("Circle.reconstructModel", (0, Circle.reconstructModel)) :+      ("Circle.hmmTrainedSupervised", (0, Circle.hmmTrainedSupervised)) :+      ("Circle.hmmTrainedUnsupervised", (0, Circle.hmmTrainedUnsupervised)) :+      ("Circle.hmmIterativelyTrained", (40, Circle.hmmIterativelyTrained)) :+      []) +++   []

@@ -0,0 +1,72 @@+module Math.HiddenMarkovModel.Utility where++import qualified Numeric.LAPACK.Matrix.Triangular as Triangular+import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape+import qualified Numeric.LAPACK.Matrix.Square as Square+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix.Triangular (Diagonal)+import Numeric.LAPACK.Vector (Vector)++import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable as StorableArray+import qualified Data.Array.Comfort.Boxed as Array+import qualified Data.Array.Comfort.Shape as Shape++import Foreign.Storable (Storable)++import qualified System.Random as Rnd++import qualified Control.Monad.Trans.State as MS+++type SquareMatrix sh = Square.Square sh++normalizeProb :: (Shape.C sh, Class.Real a) => Vector sh a -> Vector sh a+normalizeProb = snd . normalizeFactor++normalizeFactor :: (Shape.C sh, Class.Real a) => Vector sh a -> (a, Vector sh a)+normalizeFactor xs =+   let c = Vector.sum xs+   in  (c, Vector.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)+++vectorDim :: Shape.C sh => Vector sh a -> Int+vectorDim = Shape.size . StorableArray.shape++singleton :: (Class.Real a) => a -> Vector () a+singleton = Vector.constant ()+++hermitianFromList ::+   (Shape.C sh, Storable a) => sh -> [a] -> Hermitian.Hermitian sh a+hermitianFromList = Hermitian.fromList MatrixShape.RowMajor+++squareConstant ::+   (Shape.C sh, Class.Real a) => sh -> a -> SquareMatrix sh a+squareConstant = Vector.constant . MatrixShape.square MatrixShape.RowMajor++squareFromLists ::+   (Shape.C sh, Eq sh, Storable a) => sh -> [Vector sh a] -> SquareMatrix sh a+squareFromLists sh =+   Square.fromGeneral . Matrix.fromRowArray sh . Array.fromList sh++diagonal :: (Shape.C sh, Class.Real a) => Vector sh a -> Diagonal sh a+diagonal = Triangular.diagonal MatrixShape.RowMajor

@@ -0,0 +1,10 @@+module Main where++import Math.HiddenMarkovModel.Test (tests)++import qualified Test.QuickCheck as QC+++main :: IO ()+main =+   mapM_ (\(name,prop) -> putStr (name ++ ": ") >> QC.quickCheck prop) tests