hmm-lapack 0.4 → 0.4.1

raw patch · 29 files changed

+2388/−2106 lines, 29 filesdep +comfort-array-shapedep +doctest-exitcode-stdiodep +doctest-libdep ~QuickCheckdep ~basedep ~comfort-arrayPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: comfort-array-shape, doctest-exitcode-stdio, doctest-lib

Dependency ranges changed: QuickCheck, base, comfort-array, containers, deepseq, explicit-exception, fixed-length, lapack, lazy-csv, netlib-ffi, non-empty, random, semigroups, tfp, transformers, utility-ht

API changes (from Hackage documentation)

Files

• Changes.md +8/−0
• hmm-lapack.cabal +63/−21
• private/Math/HiddenMarkovModel/CSV.hs +159/−0
• private/Math/HiddenMarkovModel/Example/CirclePrivate.hs +160/−0
• private/Math/HiddenMarkovModel/Example/SineWavePrivate.hs +95/−0
• private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs +176/−0
• private/Math/HiddenMarkovModel/Normalized.hs +164/−0
• private/Math/HiddenMarkovModel/Private.hs +331/−0
• private/Math/HiddenMarkovModel/Public.hs +190/−0
• private/Math/HiddenMarkovModel/Public/Distribution.hs +548/−0
• private/Math/HiddenMarkovModel/Utility.hs +90/−0
• src/Math/HiddenMarkovModel.hs +2/−187
• src/Math/HiddenMarkovModel/CSV.hs +0/−159
• src/Math/HiddenMarkovModel/Distribution.hs +2/−539
• src/Math/HiddenMarkovModel/Example/CirclePrivate.hs +0/−123
• src/Math/HiddenMarkovModel/Example/SineWave.hs +2/−80
• src/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs +0/−157
• src/Math/HiddenMarkovModel/Normalized.hs +0/−157
• src/Math/HiddenMarkovModel/Pattern.hs +1/−1
• src/Math/HiddenMarkovModel/Private.hs +0/−331
• src/Math/HiddenMarkovModel/Test.hs +0/−259
• src/Math/HiddenMarkovModel/Utility.hs +0/−88
• test/Main.hs +9/−4
• test/Test.hs +221/−0
• test/Test/Main.hs +16/−0
• test/Test/Math/HiddenMarkovModel/Example/CirclePrivate.hs +59/−0
• test/Test/Math/HiddenMarkovModel/Example/SineWavePrivate.hs +26/−0
• test/Test/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs +48/−0
• test/Test/Math/HiddenMarkovModel/Normalized.hs +18/−0

@@ -1,3 +1,11 @@+## 0.4++* `Distribution`: Make type-classes single parameter using data families.++## 0.3++* Move to new packages `comfort-array` and `lapack`.+ ## 0.1  * `Distribution.Estimate` turned into a multi-parameter type class.

@@ -1,5 +1,6 @@+Cabal-Version:       2.2 Name:                hmm-lapack-Version:             0.4+Version:             0.4.1 Synopsis:            Hidden Markov Models using LAPACK primitives Description:   Hidden Markov Models implemented using LAPACK data types and operations.@@ -29,25 +30,24 @@   * "Math.HiddenMarkovModel.Example.Circle"   .   An alternative package without foreign calls is @hmm@.-Homepage:            http://hub.darcs.net/thielema/hmm-lapack-License:             BSD3+Homepage:            https://hub.darcs.net/thielema/hmm-lapack+License:             BSD-3-Clause 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.4+  Tag:         0.4.1   Type:        darcs-  Location:    http://hub.darcs.net/thielema/hmm-lapack+  Location:    https://hub.darcs.net/thielema/hmm-lapack  Source-Repository head   Type:        darcs-  Location:    http://hub.darcs.net/thielema/hmm-lapack+  Location:    https://hub.darcs.net/thielema/hmm-lapack  Library   Exposed-Modules:@@ -58,43 +58,85 @@     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.3 && <0.4,+    private,+    lapack >=0.4 && <0.5,     fixed-length >=0.2.1 && <0.3,     tfp >=1.0 && <1.1,     netlib-ffi >=0.1.1 && <0.2,-    comfort-array >=0.4 && <0.5,-    QuickCheck >=2.5 && <3,+    comfort-array-shape >=0.0 && <0.1,+    comfort-array >=0.5 && <0.6,     explicit-exception >=0.1.7 && <0.2,     lazy-csv >=0.5 && <0.6,-    random >=1.0 && <1.2,     transformers >= 0.2 && <0.6,     non-empty >=0.3.2 && <0.4,     semigroups >=0.17 && <1.0,     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 +Library private+  Exposed-Modules:+    Math.HiddenMarkovModel.Public+    Math.HiddenMarkovModel.Public.Distribution+    Math.HiddenMarkovModel.Example.TrafficLightPrivate+    Math.HiddenMarkovModel.Example.SineWavePrivate+    Math.HiddenMarkovModel.Example.CirclePrivate+    Math.HiddenMarkovModel.Normalized+    Math.HiddenMarkovModel.Private+    Math.HiddenMarkovModel.Utility+    Math.HiddenMarkovModel.CSV+  Build-Depends:+    lapack,+    tfp,+    netlib-ffi,+    comfort-array,+    explicit-exception,+    lazy-csv,+    random >=1.0 && <1.3,+    transformers,+    non-empty,+    semigroups,+    containers,+    utility-ht,+    deepseq,+    prelude-compat >=0.0 && <0.1,+    base+  Hs-Source-Dirs:      private+  Default-Language:    Haskell2010+  GHC-Options:         -Wall+ Test-Suite hmm-test   Type: exitcode-stdio-1.0   Build-Depends:     hmm-lapack,-    QuickCheck,+    private,+    lapack,+    comfort-array-shape,+    comfort-array,+    fixed-length,+    tfp,+    doctest-exitcode-stdio >=0.0 && <0.1,+    doctest-lib >=0.1 && <0.1.1,+    QuickCheck >=2.5 && <3,+    random,+    containers,+    non-empty,+    utility-ht,+    deepseq,     base   Main-Is: Main.hs+  Other-Modules:+    Test+    Test.Main+    Test.Math.HiddenMarkovModel.Example.TrafficLightPrivate+    Test.Math.HiddenMarkovModel.Example.SineWavePrivate+    Test.Math.HiddenMarkovModel.Example.CirclePrivate+    Test.Math.HiddenMarkovModel.Normalized   Hs-Source-Dirs:      test   Default-Language:    Haskell2010   GHC-Options:         -Wall

@@ -0,0 +1,159 @@+module Math.HiddenMarkovModel.CSV where++import Math.HiddenMarkovModel.Utility (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 (ShapeInt)+import Numeric.LAPACK.Vector (Vector)++import qualified Numeric.Netlib.Class as Class++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) => Matrix.Square 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 ShapeInt a)+parseVectorCells =+   parseVectorFields =<< getRow++-- ToDo: Maybe check row consistency already here?+parseVectorFields ::+   (Read a, Class.Real a) =>+   CSV.CSVRow -> CSVParser (Vector ShapeInt a)+parseVectorFields =+   MT.lift . fmap Vector.autoFromList . mapM parseNumberCell .+   Rev.dropWhile (null . CSV.csvFieldContent)++parseNonEmptyVectorCells ::+   (Read a, Class.Real a) =>+   CSVParser (Vector ShapeInt 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 (Matrix.Square 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 $+      Matrix.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,160 @@+module Math.HiddenMarkovModel.Example.CirclePrivate where++import qualified Math.HiddenMarkovModel.Public as HMM+import qualified Math.HiddenMarkovModel.Public.Distribution as Distr+import Math.HiddenMarkovModel.Utility+         (normalizeProb, squareFromLists, hermitianFromList)++import qualified Numeric.LAPACK.Matrix.HermitianPositiveDefinite as HermitianPD+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)+++{- $setup+>>> import qualified Math.HiddenMarkovModel as HMM+>>> import qualified Data.NonEmpty as NonEmpty+>>> import Data.Eq.HT (equating)+>>>+>>> checkTraining :: (Int, HMM) -> Bool+>>> checkTraining (maxDiff,hmm_) =+>>>    maxDiff >=+>>>    (length $ filter id $ NonEmpty.flatten $+>>>     NonEmpty.zipWith (/=)+>>>       (HMM.reveal hmm_ circle) (fmap fst circleLabeled))+-}+++data State = Q1 | Q2 | Q3 | Q4+   deriving (Eq, Ord, Enum, Bounded, Show)++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++{- |+prop> checkTraining (0, hmm)+-}+hmm :: HMM+hmm =+   HMM.Cons {+      HMM.initial = normalizeProb $ Vector.one stateSet,+      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 hermitianPD =+               HermitianPD.assurePositiveDefiniteness .+               hermitianFromList coordinateSet+             cov0 = hermitianPD [0.10, -0.09, 0.10]+             cov1 = hermitianPD [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.5+) 0++circle :: NonEmpty.T [] (Vector CoordinateSet Double)+circle = fmap snd circleLabeled++{- |+>>> take 20 $ NonEmpty.flatten revealed+[Q1,Q1,Q1,Q1,Q2,Q2,Q2,Q3,Q3,Q3,Q4,Q4,Q4,Q1,Q1,Q1,Q2,Q2,Q2,Q3]++prop> equating (take 1000 . NonEmpty.flatten) revealed $ fmap fst 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.++prop> checkTraining (0, reconstructModel)+-}+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..])+++{- |+prop> checkTraining (0, hmmTrainedSupervised)+-}+hmmTrainedSupervised :: HMM+hmmTrainedSupervised =+   HMM.finishTraining $ HMM.trainSupervised stateSet circleLabeled++{- |+prop> checkTraining (0, hmmTrainedUnsupervised)+-}+hmmTrainedUnsupervised :: HMM+hmmTrainedUnsupervised =+   HMM.finishTraining $ HMM.trainUnsupervised hmm circle++{- |+prop> checkTraining (40, hmmIterativelyTrained)+-}+hmmIterativelyTrained :: HMM+hmmIterativelyTrained =+   nest 100+      (HMM.finishTraining . flip HMM.trainUnsupervised circle)+      hmm

@@ -0,0 +1,95 @@+module Math.HiddenMarkovModel.Example.SineWavePrivate where++import qualified Math.HiddenMarkovModel.Public as HMM+import qualified Math.HiddenMarkovModel.Public.Distribution as Distr+import Math.HiddenMarkovModel.Utility (normalizeProb, squareFromLists)++import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian+import qualified Numeric.LAPACK.Matrix.Layout as Layout+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Vector (Vector, singleton)++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)+++{- $setup+>>> import qualified Data.NonEmpty as NonEmpty+-}+++data State = Rising | High | Falling | Low+   deriving (Eq, Ord, Enum, Bounded, Show)++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.one stateSet,+      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 , Hermitian.identity Layout.RowMajor ()) :+            (singleton   1 , Hermitian.identity Layout.RowMajor ()) :+            (singleton   0 , Hermitian.identity Layout.RowMajor ()) :+            (singleton (-1), Hermitian.identity Layout.RowMajor ()) :+            []+   }++stateVector :: Double -> Double -> Double -> Double -> Vector StateSet Double+stateVector x0 x1 x2 x3 = Vector.fromList stateSet [x0,x1,x2,x3]++{- |+>>> take 20 $ map fst $ NonEmpty.flatten sineWaveLabeled+[Rising,Rising,High,High,High,Falling,Falling,Falling,Low,Low,Low,Rising,Rising,Rising,Rising,High,High,High,Falling,Falling]+-}+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.5+) 0++sineWave :: NonEmpty.T [] Double+sineWave = fmap snd sineWaveLabeled++{- |+>>> take 20 $ NonEmpty.flatten revealed+[Rising,Rising,High,High,High,Falling,Falling,Falling,Low,Low,Low,Low,Rising,Rising,Rising,High,High,High,Falling,Falling]+-}+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,176 @@+{-# LANGUAGE TypeFamilies #-}+module Math.HiddenMarkovModel.Example.TrafficLightPrivate where++import qualified Math.HiddenMarkovModel.Public as HMM+import qualified Math.HiddenMarkovModel.Public.Distribution as Distr+import Math.HiddenMarkovModel.Utility (normalizeProb, squareFromLists)++import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Vector (Vector)++import qualified Data.Array.Comfort.Shape as Shape++import Text.Read.HT (maybeRead)++import Control.DeepSeq (NFData(rnf))++import qualified Data.NonEmpty as NonEmpty+import qualified Data.List.HT as ListHT+import Data.NonEmpty ((!:))+++{- $setup+>>> import qualified Data.NonEmpty as NonEmpty+>>> import Control.DeepSeq (deepseq)+>>>+>>> verifyRevelations :: HMM -> [Bool]+>>> verifyRevelations hmm_ =+>>>    map (verifyRevelation hmm_) (NonEmpty.flatten labeledSequences)+-}+++data Color = Red | Yellow | Green+   deriving (Eq, Ord, Enum, Show, Read)++instance NFData Color where+   rnf Red = ()+   rnf _ = ()++{- |+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 | 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++{- |+>>> verifyRevelations hmm+[True,True]+-}+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.discreteFromList $+            (Red,    stateVector 1 0 0 0) !:+            (Yellow, stateVector 0 1 0 1) :+            (Green,  stateVector 0 0 1 0) :+            []+   }+++{- |+>>> verifyRevelations hmmDisturbed+[True,True]+-}+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.discreteFromList $+            (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.++>>> verifyRevelations hmmTrainedSupervised+[True,True]+-}+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.++>>> verifyRevelations hmmTrainedUnsupervised+[True,True]+-}+hmmTrainedUnsupervised :: HMM+hmmTrainedUnsupervised =+   HMM.trainMany (HMM.trainUnsupervised hmm) stateSequences++{- |+Repeat unsupervised training until convergence.++prop> deepseq hmmIterativelyTrained True+-}+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)

@@ -0,0 +1,164 @@+{-# 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.Public.Distribution as Distr+import Math.HiddenMarkovModel.Private+          (T(..), Trained(..), emission,+           biscaleTransition, revealGen, sumTransitions)+import Math.HiddenMarkovModel.Utility (normalizeFactor, normalizeProb)++import qualified Numeric.LAPACK.Matrix as Matrix+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.Shape as Shape++import qualified Data.NonEmpty.Class as NonEmptyC+import qualified Data.NonEmpty as NonEmpty+import qualified Data.Foldable as Fold+import Data.Traversable (Traversable)+++{- $setup+>>> import qualified Data.NonEmpty as NonEmpty+-}+++{- |+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 typ, Shape.C sh, Eq sh, Floating prob,+    Class.Real prob, Distr.Emission typ prob ~ emission,+    Traversable f) =>+   T typ sh prob -> NonEmpty.T f emission -> prob+logLikelihood hmm = Fold.sum . fmap (log . fst) . alpha hmm++alpha ::+   (Distr.EmissionProb typ, Shape.C sh, Eq sh,+    Class.Real prob, Distr.Emission typ prob ~ emission,+    Traversable f) =>+   T typ 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 typ, Shape.C sh, Eq sh,+    Class.Real prob, Distr.Emission typ prob ~ emission,+    Traversable f, NonEmptyC.Reverse f) =>+   T typ 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.one $ StorableArray.shape $ initial hmm)++alphaBeta ::+   (Distr.EmissionProb typ, Shape.C sh, Eq sh,+    Class.Real prob, Distr.Emission typ prob ~ emission,+    Traversable f, NonEmptyC.Zip f, NonEmptyC.Reverse f) =>+   T typ 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 typ, Shape.C sh, Eq sh,+    Class.Real prob, Distr.Emission typ prob ~ emission,+    Traversable f, NonEmptyC.Zip f) =>+   T typ sh prob ->+   NonEmpty.T f emission ->+   NonEmpty.T f (prob, Vector sh prob) ->+   NonEmpty.T f (Vector sh prob) ->+   f (Matrix.Square sh prob)+xiFromAlphaBeta hmm xs calphas betas =+   let (cs,alphas) = Functor.unzip calphas+   in  NonEmptyC.zipWith4+          (\x alpha0 c1 beta1 ->+             Matrix.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 typ, Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ state,+    Distr.Emission typ prob ~ emission, Class.Real prob, Traversable f) =>+   T typ sh prob -> NonEmpty.T f emission -> NonEmpty.T f state+reveal = revealGen normalizeProb+++{- |+Variant of NonEmpty.scanr with less stack consumption.++prop> \x xs -> nonEmptyScanr (-) x xs == NonEmpty.scanr (-) x (xs::[Int])+-}+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 typ, Shape.C sh, Eq sh,+    Class.Real prob, Distr.Emission typ prob ~ emission) =>+   T typ sh prob -> NonEmpty.T [] emission -> Trained typ 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.accumulateEmissionVectors $ NonEmptyC.zip xs zetas+       }

@@ -0,0 +1,331 @@+{-# LANGUAGE TypeFamilies #-}+module Math.HiddenMarkovModel.Private where++import qualified Math.HiddenMarkovModel.Public.Distribution as Distr+import qualified Math.HiddenMarkovModel.CSV as HMMCSV+import Math.HiddenMarkovModel.Utility (diagonal)++import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix.Square as Square+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.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.Semigroup ((<>))+import Data.Traversable (Traversable, mapAccumL)+import Data.Tuple.HT (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\".+-}+data T typ sh prob =+   Cons {+      initial :: Vector sh prob,+      transition :: Matrix.Square sh prob,+      distribution :: Distr.T typ sh prob+   }+   deriving (Show)++instance+   (Distr.NFData typ, NFData sh, Shape.C sh, NFData prob, Storable prob) =>+      NFData (T typ sh prob) where+   rnf (Cons initial_ transition_ distribution_) =+      rnf (initial_, transition_, distribution_)++instance+   (Distr.Format typ, Format.FormatArray sh, Class.Real prob) =>+      Format.Format (T typ sh prob) where+   format fmt (Cons initial_ transition_ distribution_) =+      Format.format fmt (initial_, transition_, distribution_)++mapStatesShape ::+   (Distr.EmissionProb typ, Shape.C sh0, Shape.C sh1) =>+   (sh0 -> sh1) -> T typ sh0 prob -> T typ sh1 prob+mapStatesShape f hmm =+   Cons {+      initial = StorableArray.mapShape f $ initial hmm,+      transition = Square.mapSize f $ transition hmm,+      distribution = Distr.mapStatesShape f $ distribution hmm+   }+++emission ::+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob) =>+   T typ sh prob -> Distr.Emission typ prob -> Vector sh prob+emission  =  Distr.emissionProb . distribution+++forward ::+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,+    Distr.Emission typ prob ~ emission, Traversable f) =>+   T typ sh prob -> NonEmpty.T f emission -> prob+forward hmm = Vector.sum . NonEmpty.last . alpha hmm++alpha ::+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,+    Distr.Emission typ prob ~ emission, Traversable f) =>+   T typ 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 ::+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,+    Distr.Emission typ prob ~ emission, Traversable f) =>+   T typ sh prob -> NonEmpty.T f emission -> prob+backward hmm (NonEmpty.Cons x xs) =+   Vector.dot (initial hmm) $+   Vector.mul (emission hmm x) $+   NonEmpty.head $ beta hmm xs++beta ::+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,+    Distr.Emission typ prob ~ emission, Traversable f) =>+   T typ 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.one $ StorableArray.shape $ initial hmm)+++alphaBeta ::+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,+    Distr.Emission typ prob ~ emission, Traversable f) =>+   T typ 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 ::+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob) =>+   T typ sh prob -> Distr.Emission typ prob ->+   Vector sh prob -> Vector sh prob -> Matrix.Square sh prob+biscaleTransition hmm x alpha0 beta1 =+   (diagonal (Vector.mul (emission hmm x) beta1)+    #*##+    transition hmm)+   ##*#+   diagonal alpha0++xiFromAlphaBeta ::+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,+    Distr.Emission typ prob ~ emission) =>+   T typ sh prob -> prob ->+   NonEmpty.T [] emission ->+   NonEmpty.T [] (Vector sh prob) ->+   NonEmpty.T [] (Vector sh prob) ->+   [Matrix.Square sh prob]+xiFromAlphaBeta hmm recipLikelihood xs alphas betas =+   zipWith3+      (\x alpha0 beta1 ->+         Matrix.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) =>+   [Matrix.Square 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 ::+   (Distr.EmissionProb typ, Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ state,+    Distr.Emission typ prob ~ emission, Class.Real prob, Traversable f) =>+   T typ sh prob -> NonEmpty.T f emission -> NonEmpty.T f state+reveal = revealGen id++revealGen ::+   (Distr.EmissionProb typ, Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ state,+    Distr.Emission typ prob ~ emission, Class.Real prob, Traversable f) =>+   (Vector (Shape.Deferred sh) prob -> Vector (Shape.Deferred sh) prob) ->+   T typ sh prob -> NonEmpty.T f emission -> NonEmpty.T f state+revealGen normalize hmm =+   fmap (Shape.revealIndex (StorableArray.shape $ initial hmm)) .+   revealStorable normalize (mapStatesShape Shape.Deferred hmm)++revealStorable ::+   (Distr.EmissionProb typ, Shape.InvIndexed sh, Eq sh,+    Shape.Index sh ~ state, Storable state,+    Distr.Emission typ prob ~ emission, Class.Real prob, Traversable f) =>+   (Vector sh prob -> Vector sh prob) ->+   T typ sh prob -> NonEmpty.T f emission -> NonEmpty.T f state+revealStorable normalize hmm (NonEmpty.Cons x xs) =+   uncurry (NonEmpty.scanr (StorableArray.!)) $+   mapFst (fst . Vector.argAbsMaximum) $+   mapAccumL+      (\alphai xi ->+         swap $ mapSnd (Vector.mul (emission hmm xi)) $+         matrixMaxMul (transition hmm) $ normalize alphai)+      (Vector.mul (emission hmm x) (initial hmm)) xs++matrixMaxMul ::+   (Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ ix, Storable ix,+    Class.Real a) =>+   Matrix.Square sh a -> Vector sh a ->+   (Vector sh ix, Vector sh a)+matrixMaxMul m v = Matrix.rowArgAbsMaximums $ Matrix.scaleColumns v 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 typ sh prob =+   Trained {+      trainedInitial :: Vector sh prob,+      trainedTransition :: Matrix.Square sh prob,+      trainedDistribution :: Distr.Trained typ sh prob+   }+   deriving (Show)++instance+   (Distr.NFData typ, NFData sh, Shape.C sh, NFData prob, Storable prob) =>+      NFData (Trained typ sh prob) where+   rnf hmm =+      rnf (trainedInitial hmm, trainedTransition hmm, trainedDistribution hmm)+++sumTransitions ::+   (Shape.C sh, Eq sh, Class.Real e) =>+   T typ sh e -> [Matrix.Square sh e] -> Matrix.Square sh e+sumTransitions hmm =+   List.foldl' Matrix.add $+   Matrix.zero $ ArrMatrix.shape $ transition hmm++{- |+Baum-Welch algorithm+-}+trainUnsupervised ::+   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob,+    Distr.Emission typ prob ~ emission) =>+   T typ sh prob -> NonEmpty.T [] emission -> Trained typ 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.accumulateEmissionVectors $ NonEmptyC.zip xs zetas+       }+++mergeTrained ::+   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>+   Trained typ sh prob -> Trained typ sh prob -> Trained typ sh prob+mergeTrained hmm0 hmm1 =+   Trained {+      trainedInitial = Vector.add (trainedInitial hmm0) (trainedInitial hmm1),+      trainedTransition =+         Matrix.add (trainedTransition hmm0) (trainedTransition hmm1),+      trainedDistribution =+         trainedDistribution hmm0 <> trainedDistribution hmm1+   }++instance+   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>+      Sg.Semigroup (Trained typ sh prob) where+   (<>) = mergeTrained+++toCells ::+   (Distr.ToCSV typ, Shape.Indexed sh, Class.Real prob, Show prob) =>+   T typ sh prob -> [[String]]+toCells hmm =+   (HMMCSV.cellsFromVector $ initial hmm) :+   (HMMCSV.cellsFromSquare $ transition hmm) +++   [] :+   (Distr.toCells $ distribution hmm)++parseCSV ::+   (Distr.FromCSV typ, Shape.C stateSh, Eq stateSh,+    Class.Real prob, Read prob) =>+   (Int -> stateSh) -> HMMCSV.CSVParser (T typ 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,190 @@+{-# LANGUAGE TypeFamilies #-}+module Math.HiddenMarkovModel.Public (+   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.Public.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+          (squareConstant, distance, matrixDistance,+           randomItemProp, normalizeProb, attachOnes)++import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix ((#!))++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 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.Discrete symbol) sh prob+type Discrete symbol sh prob = T (Distr.Discrete symbol) sh prob++type GaussianTrained emiSh stateSh a =+         Trained (Distr.Gaussian emiSh) stateSh a+type Gaussian emiSh stateSh a = T (Distr.Gaussian emiSh) 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 typ, Shape.C sh, Class.Real prob) =>+   Distr.T typ sh prob -> T typ 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 ::+   (Distr.EmissionProb typ, Shape.Indexed sh, Shape.Index sh ~ state,+    Class.Real prob, Distr.Emission typ prob ~ emission, Traversable f) =>+   T typ sh prob -> f (state, emission) -> f prob+probabilitySequence hmm =+   snd+   .+   mapAccumL+      (\index (s, e) ->+         ((transition hmm #!) . flip (,) s,+          index s * Distr.emissionStateProb (distribution hmm) e s))+      (initial hmm StorableArray.!)++generate ::+   (Distr.Generate typ, Shape.Indexed sh, Class.Real prob,+    Rnd.RandomGen g, Rnd.Random prob, Distr.Emission typ prob ~ emission) =>+   T typ sh prob -> g -> [emission]+generate hmm = map snd . generateLabeled hmm++generateLabeled ::+   (Distr.Generate typ, Shape.Indexed sh, Shape.Index sh ~ state,+    Rnd.RandomGen g, Rnd.Random prob,+    Class.Real prob, Distr.Emission typ prob ~ emission) =>+   T typ 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.Estimate typ, Shape.Indexed sh, Shape.Index sh ~ state,+    Class.Real prob, Distr.Emission typ prob ~ emission) =>+   sh -> NonEmpty.T [] (state, emission) -> Trained typ sh prob+trainSupervised sh xs =+   let getState (s, _x) = s+   in  Trained {+          trainedInitial = Vector.unit sh $ getState $ NonEmpty.head xs,+          trainedTransition =+             Matrix.transpose $ ArrMatrix.fromVector $+             StorableArray.accumulate (+)+                (ArrMatrix.toVector $ squareConstant sh 0) $+             attachOnes $ NonEmpty.mapAdjacent (,) $ fmap getState xs,+          trainedDistribution = Distr.accumulateEmissions sh xs+       }++finishTraining ::+   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>+   Trained typ sh prob -> T typ 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) => Matrix.Square sh a -> Matrix.Square sh a+normalizeProbColumns m =+   Matrix.scaleColumns (StorableArray.map recip (Matrix.columnSums m)) m++trainMany ::+   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob, Foldable f) =>+   (trainingData -> Trained typ sh prob) ->+   NonEmpty.T f trainingData -> T typ 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.C sh, Eq sh, Class.Real prob) =>+   T typ sh prob -> T typ sh prob -> prob+deviation hmm0 hmm1 =+   distance (initial hmm0) (initial hmm1)+   `max`+   matrixDistance (transition hmm0) (transition hmm1)+++toCSV ::+   (Distr.ToCSV typ, Shape.Indexed sh, Class.Real prob, Show prob) =>+   T typ sh prob -> String+toCSV hmm =+   CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $ toCells hmm++fromCSV ::+   (Distr.FromCSV typ, Shape.Indexed sh, Eq sh, Class.Real prob, Read prob) =>+   (Int -> sh) -> String -> ME.Exceptional String (T typ sh prob)+fromCSV makeShape =+   MS.evalStateT (parseCSV makeShape) . map HMMCSV.fixShortRow . CSV.parseCSV

@@ -0,0 +1,548 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE EmptyDataDecls #-}+module Math.HiddenMarkovModel.Public.Distribution (+   T(..), Trained(..), Emission,+   Show(..), NFData(..), Format(..),+   Info(..), Generate(..), EmissionProb(..),+   Estimate(..), accumulateEmissionVectors,++   Discrete, discreteFromList,+   Gaussian, gaussian, gaussianTrained,++   ToCSV(..), FromCSV(..), HMMCSV.CSVParser, CSVSymbol(..),+   ) where++import qualified Math.HiddenMarkovModel.CSV as HMMCSV+import Math.HiddenMarkovModel.Utility (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.Layout as Layout+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import qualified Numeric.LAPACK.Format as Format+import qualified Numeric.LAPACK.Output as Output+import Numeric.LAPACK.Matrix ((-*#), (-/#), (#/\), (|*-), (#!))+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Format (FormatArray)+import Numeric.LAPACK.Output (Output)++import qualified Type.Data.Bool as TBool++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 Data.Array.Comfort.Shape ((::+)((::+)))++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 qualified Control.DeepSeq as DeepSeq+import Control.Monad (liftM2)+import Control.Applicative (liftA2)++import qualified Data.NonEmpty.Map as NonEmptyMap+import qualified Data.NonEmpty as NonEmpty+import qualified Data.Semigroup as Sg+import qualified Data.Map as Map+import qualified Data.List.HT as ListHT+import qualified Data.List as List+import Data.Functor.Identity (Identity(Identity), runIdentity)+import Data.Tuple.HT (snd3)+import Data.Set (Set)+import Data.Maybe (listToMaybe)++import qualified Prelude as P+import Prelude2010 hiding (Show, showsPrec)++++data family T typ sh prob+data family Trained typ sh prob++type family Emission typ prob+++class Show typ where+   showsPrec ::+      (Shape.C sh, P.Show sh, P.Show prob, Storable prob) =>+      Int -> T typ sh prob -> ShowS+   showsPrecTrained ::+      (Shape.C sh, P.Show sh, P.Show prob, Storable prob) =>+      Int -> Trained typ sh prob -> ShowS++instance+   (Show typ, Shape.C sh, P.Show sh, P.Show prob, Storable prob) =>+      P.Show (T typ sh prob) where+   showsPrec = showsPrec++instance+   (Show typ, Shape.C sh, P.Show sh, P.Show prob, Storable prob) =>+      P.Show (Trained typ sh prob) where+   showsPrec = showsPrecTrained+++class NFData typ where+   rnf ::+      (DeepSeq.NFData sh, DeepSeq.NFData prob, Shape.C sh) =>+      T typ sh prob -> ()+   rnfTrained ::+      (DeepSeq.NFData sh, DeepSeq.NFData prob, Shape.C sh) =>+      Trained typ sh prob -> ()++instance+   (NFData typ, DeepSeq.NFData sh, DeepSeq.NFData prob, Shape.C sh) =>+      DeepSeq.NFData (T typ sh prob) where+   rnf = rnf++instance+   (NFData typ, DeepSeq.NFData sh, DeepSeq.NFData prob, Shape.C sh) =>+      DeepSeq.NFData (Trained typ sh prob) where+   rnf = rnfTrained+++class Format typ where+   format ::+      (Shape.C sh, Output out, Class.Real prob) =>+      String -> T typ sh prob -> out++instance+   (Format typ, Shape.C sh, Class.Real prob) =>+      Format.Format (T typ sh prob) where+   format = format++++class Info typ where+   statesShape :: (Shape.C sh) => T typ sh prob -> sh+   statesShapeTrained :: (Shape.C sh) => Trained typ sh prob -> sh++class Generate typ where+   generate ::+      (Shape.Indexed sh, Class.Real prob, Rnd.Random prob, Rnd.RandomGen g) =>+      T typ sh prob -> Shape.Index sh -> MS.State g (Emission typ prob)++class EmissionProb typ where+   mapStatesShape ::+      (Shape.C sh0, Shape.C sh1) =>+      (sh0 -> sh1) -> T typ sh0 prob -> T typ sh1 prob+   {-+   This function could be implemented generically in terms of emissionStateProb+   but that would require an Info constraint.+   -}+   emissionProb ::+      (Shape.C sh, Class.Real prob) =>+      T typ sh prob -> Emission typ prob -> Vector sh prob+   emissionStateProb ::+      (Shape.Indexed sh, Class.Real prob) =>+      T typ sh prob -> Emission typ prob -> Shape.Index sh -> prob+   emissionStateProb distr e s = emissionProb distr e StorableArray.! s++class (EmissionProb typ) => Estimate typ where+   accumulateEmissions ::+      (Shape.Indexed sh, Class.Real prob, Shape.Index sh ~ state) =>+      sh -> NonEmpty.T [] (state, Emission typ prob) -> Trained typ sh prob+   trainVector ::+      (Shape.C sh, Eq sh, Class.Real prob) =>+      Emission typ prob -> Vector sh prob -> Trained typ sh prob+   combine ::+      (Shape.C sh, Eq sh, Class.Real prob) =>+      Trained typ sh prob -> Trained typ sh prob -> Trained typ sh prob+   normalize ::+      (Shape.C sh, Eq sh, Class.Real prob) =>+      Trained typ sh prob -> T typ sh prob++accumulateEmissionVectors ::+   (Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>+   NonEmpty.T [] (Emission typ prob, Vector sh prob) -> Trained typ sh prob+accumulateEmissionVectors = NonEmpty.foldl1Map combine (uncurry trainVector)++instance+   (Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>+      Sg.Semigroup (Trained typ sh prob) where+   (<>) = combine+++data Discrete symbol++newtype instance T (Discrete symbol) sh prob =+      Discrete (Matrix.General (Set symbol) sh prob)++newtype instance Trained (Discrete symbol) sh prob =+      DiscreteTrained (NonEmptyMap.T symbol (Vector sh prob))++type instance Emission (Discrete symbol) prob = symbol+++instance (P.Show symbol, Ord symbol) => Show (Discrete symbol) where+   showsPrec prec (Discrete m) = P.showsPrec prec m+   showsPrecTrained prec (DiscreteTrained m) = P.showsPrec prec m++instance (DeepSeq.NFData symbol) => NFData (Discrete symbol) where+   rnf (Discrete m) = DeepSeq.rnf m+   rnfTrained (DiscreteTrained m) = DeepSeq.rnf m++instance (P.Show symbol, Ord symbol) => Format (Discrete symbol) where+   format fmt (Discrete m) =+      Output.formatAligned $+      map (\(sym,v) ->+            map (Identity . Output.text) $+            (show sym ++ ":") : map (printFmt fmt) (Vector.toList v)) $+      Array.toAssociations $ Matrix.toRowArray m++-- cf. Data.Bifunctor.Flip+newtype Flip f b a = Flip {getFlip :: f a b}++printFmt :: (Class.Real a) => String -> a -> String+printFmt fmt =+   getFlip $ Class.switchReal (Flip $ printf fmt) (Flip $ printf fmt)++instance (Ord symbol) => Info (Discrete symbol) where+   statesShape (Discrete m) = Matrix.width m+   statesShapeTrained (DiscreteTrained m) = discreteStateShape m++instance (Ord symbol) => Generate (Discrete symbol) where+   generate (Discrete m) =+      randomItemProp . StorableArray.toAssociations . Matrix.takeColumn m++instance (Ord symbol) => EmissionProb (Discrete symbol) where+   mapStatesShape f (Discrete m) = Discrete $ Matrix.mapWidth f m+   emissionProb (Discrete m) = Matrix.takeRow m+   emissionStateProb (Discrete m) x s = m #! (x,s)++instance (Ord symbol) => Estimate (Discrete symbol) where+   accumulateEmissions sh =+      DiscreteTrained .+      NonEmptyMap.map+         (StorableArray.reshape sh .+          StorableArray.fromAssociations 0 (Shape.Deferred sh) .+          Map.toList) .+      NonEmptyMap.fromListWith (Map.unionWith (+)) .+      fmap (\(state,sym) -> (sym, Map.singleton (Shape.deferIndex sh state) 1))+   trainVector sym = DiscreteTrained . NonEmptyMap.singleton sym+   combine (DiscreteTrained distr0) (DiscreteTrained distr1) =+      DiscreteTrained $ NonEmptyMap.unionWith Vector.add distr0 distr1+   normalize (DiscreteTrained distr) =+      Discrete $ normalizeProbColumns $ discreteFromMap distr++normalizeProbColumns ::+   (Shape.C height, Shape.C width, Eq width, Class.Real a) =>+   Matrix.General height width a -> Matrix.General height width a+normalizeProbColumns m = m #/\ Matrix.columnSums m++discreteStateShape ::+   (Shape.C sh) => NonEmptyMap.T symbol (Vector sh prob) -> sh+discreteStateShape =+   StorableArray.shape . snd . fst . NonEmptyMap.minViewWithKey++discreteFromMap ::+   (Ord symbol, Shape.C sh, Eq sh, Class.Real prob) =>+   NonEmptyMap.T symbol (Vector sh prob) -> Matrix.General (Set symbol) sh prob+discreteFromMap m =+   Matrix.fromRowArray (discreteStateShape m) $+   Array.fromMap $ NonEmptyMap.flatten m++discreteFromList ::+   (Ord symbol, Shape.C sh, Eq sh, Class.Real prob) =>+   NonEmpty.T [] (symbol, Vector sh prob) -> T (Discrete symbol) sh prob+discreteFromList = Discrete . discreteFromMap . NonEmptyMap.fromList++++data Gaussian emiSh++newtype instance T (Gaussian emiSh) stateSh a =+   Gaussian (Array stateSh (a, Vector emiSh a, Triangular.Upper emiSh a))++newtype instance Trained (Gaussian emiSh) stateSh a =+   GaussianTrained+      (StorableArray.Array (stateSh, Layout.Hermitian (()::+emiSh)) a)++type instance Emission (Gaussian emiSh) a = Vector emiSh a+++instance (Shape.C emiSh, P.Show emiSh) => Show (Gaussian emiSh) where+   showsPrec prec (Gaussian m) = P.showsPrec prec m+   showsPrecTrained prec (GaussianTrained m) = P.showsPrec prec m++instance (DeepSeq.NFData emiSh) => NFData (Gaussian emiSh) where+   rnf (Gaussian params) = DeepSeq.rnf params+   rnfTrained (GaussianTrained params) = DeepSeq.rnf params+++instance (FormatArray emiSh) => Format (Gaussian emiSh) where+   format = runFormatGaussian $ Class.switchReal formatGaussian formatGaussian++newtype FormatGaussian out emiSh stateSh a =+   FormatGaussian+      {runFormatGaussian :: String -> T (Gaussian emiSh) stateSh a -> out}++formatGaussian ::+   (FormatArray emiSh, Shape.C stateSh,+    Class.Real a, Format.Format a, Output out) =>+   FormatGaussian out emiSh stateSh a+formatGaussian =+   FormatGaussian $ \fmt (Gaussian params) ->+      Format.format fmt $ Array.toList params+++instance Info (Gaussian emiSh) where+   statesShape (Gaussian params) = Array.shape params+   statesShapeTrained (GaussianTrained params) =+      fst $ StorableArray.shape params++instance (Shape.C emiSh, Eq emiSh) => Generate (Gaussian emiSh) where+   generate (Gaussian allParams) state = do+      let (_c, center, covarianceChol) = 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) => EmissionProb (Gaussian emiSh) where+   mapStatesShape f (Gaussian m) = Gaussian $ Array.mapShape f m+   emissionProb (Gaussian allParams) x =+      StorableArray.fromBoxed $ fmap (gaussianEmissionProb x) allParams+   emissionStateProb (Gaussian allParams) x s =+      gaussianEmissionProb x $ allParams ! s++gaussianEmissionProb ::+   (Shape.C emiSh, Eq emiSh, Class.Real a) =>+   Vector emiSh a -> (a, Vector emiSh a, Triangular.Upper emiSh a) -> a+gaussianEmissionProb x (c, center, covarianceChol) =+   c * expSquared (Vector.sub x center -/# covarianceChol)++expSquared :: (Shape.C sh, Class.Real a) => Vector sh a -> a+expSquared =+   getNorm $ Class.switchReal (Norm expSquaredAux) (Norm expSquaredAux)++newtype Norm f a = Norm {getNorm :: f a -> a}++expSquaredAux ::+   (Shape.C sh, Class.Floating a, Vector.RealOf a ~ ar, Class.Real ar) =>+   Vector sh a -> ar+expSquaredAux x = exp ((-1/2) * Vector.norm2Squared x)+++instance (Shape.C emiSh, Eq emiSh) => Estimate (Gaussian emiSh) where+   accumulateEmissions sh xs =+      let emiSh = StorableArray.shape $ snd $ NonEmpty.head xs+          hermSh = Layout.hermitian Layout.RowMajor (()::+emiSh)+      in GaussianTrained $+         Matrix.toRowMajor . Matrix.fromRowArray hermSh . Array.reshape sh .+         Array.accumulate Vector.add+            (Array.replicate (Shape.Deferred sh) (Vector.zero hermSh)) .+         map (\(state,v) -> (Shape.deferIndex sh state, extendedHermitian v)) .+         NonEmpty.flatten+            $ xs+   trainVector xs probs =+      GaussianTrained $ Matrix.toRowMajor $ probs |*- extendedHermitian xs+   combine (GaussianTrained m0) (GaussianTrained m1) =+      GaussianTrained $ Vector.add m0 m1+   {-+     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 m) =+      let params (weight, centerSum, covarianceSum) =+             let c = recip (weight#!((),()))+                 center = Vector.scale c $ Matrix.flattenRow centerSum+             in  (center,+                  HermitianPD.assurePositiveDefiniteness $+                  Matrix.sub+                     (Matrix.scaleRealReal c covarianceSum)+                     (Hermitian.relaxIndefinite $+                      Hermitian.outer Layout.RowMajor center))+      in Gaussian $+         fmap (gaussianParameters . params .+               Hermitian.split . ArrMatrix.fromVector) $+         Matrix.toRowArray $ Matrix.fromRowMajor m++extendedHermitian ::+   (Shape.C emiSh, Class.Floating a) =>+   StorableArray.Array emiSh a ->+   StorableArray.Array (Layout.Hermitian (()::+emiSh)) a+extendedHermitian =+   ArrMatrix.toVector .+   Hermitian.outer Layout.RowMajor . Vector.append (Vector.one ())++{- |+input array must be non-empty+-}+gaussianTrained ::+   (TBool.C zero, Shape.C emiSh, Eq emiSh, Shape.C stateSh, Class.Real prob) =>+   Array stateSh+      (prob, Vector emiSh prob,+       Matrix.FlexHermitian TBool.False zero TBool.True emiSh prob) ->+   Trained (Gaussian emiSh) stateSh prob+gaussianTrained =+   GaussianTrained . Matrix.toRowMajor .+   matrixFromRowArray "HMM.Distribution.gaussianTrained" .+   fmap+      (\(weight, center, covariance) ->+         ArrMatrix.toVector $+         Hermitian.stack+            (Hermitian.fromList Layout.RowMajor () [weight])+            (Matrix.singleRow Layout.RowMajor center)+            (Hermitian.relaxIndefinite covariance))++matrixFromRowArray ::+   (Shape.C width, Eq width, Shape.C height, Class.Real a) =>+   String ->+   Array height (StorableArray.Array width a) ->+   Matrix.General height width a+matrixFromRowArray name xs =+   case Array.toList xs of+      [] -> error $ name ++ ": empty array"+      x:_ -> Matrix.fromRowArray (StorableArray.shape x) xs++gaussian ::+   (Shape.C emiSh, Shape.C stateSh, Class.Real prob) =>+   Array stateSh (Vector emiSh prob, Matrix.HermitianPosDef emiSh prob) ->+   T (Gaussian emiSh) stateSh prob+gaussian = Gaussian . fmap gaussianParameters++gaussianParameters ::+   (Shape.C emiSh, Class.Real prob) =>+   (Vector emiSh prob, Matrix.HermitianPosDef emiSh prob) ->+   (prob, Vector emiSh prob, Triangular.Upper emiSh prob)+gaussianParameters (center, covariance) =+   gaussianFromCholesky center $ HermitianPD.decompose covariance++gaussianFromCholesky ::+   (Shape.C emiSh, Class.Real prob) =>+   Vector emiSh prob -> Triangular.Upper emiSh prob ->+   (prob, Vector emiSh prob, Triangular.Upper emiSh prob)+gaussianFromCholesky center covarianceChol =+   let covarianceSqrtDet =+         Vector.product $ Triangular.takeDiagonal covarianceChol+   in  (recip (sqrt2pi ^ vectorDim center * covarianceSqrtDet),+        center, covarianceChol)++sqrt2pi :: (Class.Real a) => a+sqrt2pi = runIdentity $ Class.switchReal sqrt2piAux sqrt2piAux++sqrt2piAux :: (Floating a) => Identity a+sqrt2piAux = Identity $ sqrt (2*pi)+++class ToCSV typ where+   toCells ::+      (Shape.C sh, Class.Real prob, P.Show prob) =>+      T typ sh prob -> [[String]]++class FromCSV typ where+   parseCells ::+      (Shape.C sh, Eq sh, Class.Real prob, Read prob) =>+      sh -> HMMCSV.CSVParser (T typ sh prob)++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 (CSVSymbol symbol) => ToCSV (Discrete symbol) where+   toCells (Discrete m) =+      map+         (\(symbol, probs) ->+            cellFromSymbol symbol : HMMCSV.cellsFromVector probs) $+      Array.toAssociations $ Matrix.toRowArray m++instance (CSVSymbol symbol) => FromCSV (Discrete symbol) where+   parseCells n =+      let p = parseSymbolProb n+      in fmap discreteFromList $+         liftA2 NonEmpty.Cons (HMMCSV.getRow >>= p) (HMMCSV.manyRowsUntilEnd p)++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) => ToCSV (Gaussian emiSh) where+   toCells (Gaussian params) =+      List.intercalate [[]] $+      map+         (\(_, center, covarianceChol) ->+            HMMCSV.cellsFromVector center :+            HMMCSV.cellsFromSquare (Triangular.toSquare covarianceChol)) $+      Array.toList params++instance (emiSh ~ Matrix.ShapeInt) => FromCSV (Gaussian emiSh) 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 . snd3) gs+      HMMCSV.assert (ListHT.allEqual sizes) $+         printf "dimensions of emissions mismatch: %s" (show sizes)+      return $ Gaussian $ Array.fromList sh gs++parseSingleGaussian ::+   (emiSh ~ Matrix.ShapeInt, Class.Real prob, Eq prob, Read prob) =>+   HMMCSV.CSVParser (prob, Vector emiSh prob, Triangular.Upper emiSh 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) =>+   Matrix.Square sh a -> Triangular.Upper sh a -> Bool+isUpperTriang m mt =+   Vector.toList (ArrMatrix.toVector m)+   ==+   Vector.toList (ArrMatrix.toVector (Triangular.toSquare mt))

@@ -0,0 +1,90 @@+module Math.HiddenMarkovModel.Utility where++import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian+import qualified Numeric.LAPACK.Matrix.Layout as Layout+import qualified Numeric.LAPACK.Matrix.Extent as Extent+import qualified Numeric.LAPACK.Matrix.Square as Square+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix.Array (ArrayMatrix)+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+++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, 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+++hermitianFromList ::+   (Shape.C sh, Class.Floating a) => sh -> [a] -> Hermitian.Hermitian sh a+hermitianFromList = Hermitian.fromList Layout.RowMajor+++squareConstant ::+   (Shape.C sh, Class.Real a) => sh -> a -> Matrix.Square sh a+squareConstant =+   (ArrMatrix.fromVector .) .+   Vector.constant . Layout.square Layout.RowMajor++squareFromLists ::+   (Shape.C sh, Eq sh, Storable a) => sh -> [Vector sh a] -> Matrix.Square sh a+squareFromLists sh =+   Square.fromFull . Matrix.fromRowArray sh . Array.fromList sh++diagonal :: (Shape.C sh, Class.Real a) => Vector sh a -> Matrix.Diagonal sh a+diagonal = Matrix.diagonal Layout.RowMajor+++newtype Distance f a = Distance {getDistance :: f a -> f a -> a}++distance ::+   (Shape.C sh, Eq sh, Class.Real a) =>+   Vector sh a -> Vector sh a -> a+distance =+   getDistance $+   Class.switchReal+      (Distance $ (Vector.normInf .) . Vector.sub)+      (Distance $ (Vector.normInf .) . Vector.sub)++matrixDistance ::+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>+   (Shape.C height, Shape.C width, Eq height, Eq width, Class.Real a) =>+   ArrayMatrix pack prop lower upper meas vert horiz height width a ->+   ArrayMatrix pack prop lower upper meas vert horiz height width a ->+   a+matrixDistance a b = distance (ArrMatrix.unwrap a) (ArrMatrix.unwrap b)

@@ -1,190 +1,5 @@-{-# 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,+   module Math.HiddenMarkovModel.Public    ) 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-          (squareConstant, distance, matrixDistance,-           randomItemProp, normalizeProb, attachOnes)--import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix-import qualified Numeric.LAPACK.Matrix as Matrix-import qualified Numeric.LAPACK.Vector as Vector-import Numeric.LAPACK.Matrix ((#!))--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 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.Discrete symbol) sh prob-type Discrete symbol sh prob = T (Distr.Discrete symbol) sh prob--type GaussianTrained emiSh stateSh a =-         Trained (Distr.Gaussian emiSh) stateSh a-type Gaussian emiSh stateSh a = T (Distr.Gaussian emiSh) 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 typ, Shape.C sh, Class.Real prob) =>-   Distr.T typ sh prob -> T typ 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 ::-   (Distr.EmissionProb typ, Shape.Indexed sh, Shape.Index sh ~ state,-    Class.Real prob, Distr.Emission typ prob ~ emission, Traversable f) =>-   T typ sh prob -> f (state, emission) -> f prob-probabilitySequence hmm =-   snd-   .-   mapAccumL-      (\index (s, e) ->-         ((transition hmm #!) . flip (,) s,-          index s * Distr.emissionStateProb (distribution hmm) e s))-      (initial hmm StorableArray.!)--generate ::-   (Distr.Generate typ, Shape.Indexed sh, Class.Real prob,-    Rnd.RandomGen g, Rnd.Random prob, Distr.Emission typ prob ~ emission) =>-   T typ sh prob -> g -> [emission]-generate hmm = map snd . generateLabeled hmm--generateLabeled ::-   (Distr.Generate typ, Shape.Indexed sh, Shape.Index sh ~ state,-    Rnd.RandomGen g, Rnd.Random prob,-    Class.Real prob, Distr.Emission typ prob ~ emission) =>-   T typ 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.Estimate typ, Shape.Indexed sh, Shape.Index sh ~ state,-    Class.Real prob, Distr.Emission typ prob ~ emission) =>-   sh -> NonEmpty.T [] (state, emission) -> Trained typ sh prob-trainSupervised sh xs =-   let getState (s, _x) = s-   in  Trained {-          trainedInitial = Vector.unit sh $ getState $ NonEmpty.head xs,-          trainedTransition =-             Matrix.transpose $ ArrMatrix.fromVector $-             StorableArray.accumulate (+)-                (ArrMatrix.toVector $ squareConstant sh 0) $-             attachOnes $ NonEmpty.mapAdjacent (,) $ fmap getState xs,-          trainedDistribution = Distr.accumulateEmissions sh xs-       }--finishTraining ::-   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>-   Trained typ sh prob -> T typ 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) => Matrix.Square sh a -> Matrix.Square sh a-normalizeProbColumns m =-   Matrix.scaleColumns (StorableArray.map recip (Matrix.columnSums m)) m--trainMany ::-   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob, Foldable f) =>-   (trainingData -> Trained typ sh prob) ->-   NonEmpty.T f trainingData -> T typ 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.C sh, Eq sh, Class.Real prob) =>-   T typ sh prob -> T typ sh prob -> prob-deviation hmm0 hmm1 =-   distance (initial hmm0) (initial hmm1)-   `max`-   matrixDistance (transition hmm0) (transition hmm1)---toCSV ::-   (Distr.ToCSV typ, Shape.Indexed sh, Class.Real prob, Show prob) =>-   T typ sh prob -> String-toCSV hmm =-   CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $ toCells hmm--fromCSV ::-   (Distr.FromCSV typ, Shape.Indexed sh, Eq sh, Class.Real prob, Read prob) =>-   (Int -> sh) -> String -> ME.Exceptional String (T typ sh prob)-fromCSV makeShape =-   MS.evalStateT (parseCSV makeShape) . map HMMCSV.fixShortRow . CSV.parseCSV+import Math.HiddenMarkovModel.Public

@@ -1,159 +0,0 @@-module Math.HiddenMarkovModel.CSV where--import Math.HiddenMarkovModel.Utility (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 (ShapeInt)-import Numeric.LAPACK.Vector (Vector)--import qualified Numeric.Netlib.Class as Class--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) => Matrix.Square 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 ShapeInt a)-parseVectorCells =-   parseVectorFields =<< getRow---- ToDo: Maybe check row consistency already here?-parseVectorFields ::-   (Read a, Class.Real a) =>-   CSV.CSVRow -> CSVParser (Vector ShapeInt a)-parseVectorFields =-   MT.lift . fmap Vector.autoFromList . mapM parseNumberCell .-   Rev.dropWhile (null . CSV.csvFieldContent)--parseNonEmptyVectorCells ::-   (Read a, Class.Real a) =>-   CSVParser (Vector ShapeInt 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 (Matrix.Square 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 $-      Matrix.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)

@@ -1,542 +1,5 @@-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE EmptyDataDecls #-} module Math.HiddenMarkovModel.Distribution (-   T(..), Trained(..), Emission,-   Show(..), NFData(..), Format(..),-   Info(..), Generate(..), EmissionProb(..),-   Estimate(..), accumulateEmissionVectors,--   Discrete, discreteFromList,-   Gaussian, gaussian, gaussianTrained,--   ToCSV(..), FromCSV(..), HMMCSV.CSVParser, CSVSymbol(..),+   module Math.HiddenMarkovModel.Public.Distribution    ) where -import qualified Math.HiddenMarkovModel.CSV as HMMCSV-import Math.HiddenMarkovModel.Utility (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.Array as ArrMatrix-import qualified Numeric.LAPACK.Matrix as Matrix-import qualified Numeric.LAPACK.Vector as Vector-import qualified Numeric.LAPACK.Format as Format-import qualified Numeric.LAPACK.Output as Output-import Numeric.LAPACK.Matrix ((-*#), (-/#), (#/\), (|*-), (#!))-import Numeric.LAPACK.Vector (Vector)-import Numeric.LAPACK.Format (FormatArray)-import Numeric.LAPACK.Output (Output)--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 Data.Array.Comfort.Shape ((:+:)((:+:)))--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 qualified Control.DeepSeq as DeepSeq-import Control.Monad (liftM2)-import Control.Applicative (liftA2)--import qualified Data.NonEmpty.Map as NonEmptyMap-import qualified Data.NonEmpty as NonEmpty-import qualified Data.Semigroup as Sg-import qualified Data.Map as Map-import qualified Data.List.HT as ListHT-import qualified Data.List as List-import Data.Functor.Identity (Identity(Identity), runIdentity)-import Data.Tuple.HT (snd3)-import Data.Set (Set)-import Data.Maybe (listToMaybe)--import qualified Prelude as P-import Prelude2010 hiding (Show, showsPrec)----data family T typ sh prob-data family Trained typ sh prob--type family Emission typ prob---class Show typ where-   showsPrec ::-      (Shape.C sh, P.Show sh, P.Show prob, Storable prob) =>-      Int -> T typ sh prob -> ShowS-   showsPrecTrained ::-      (Shape.C sh, P.Show sh, P.Show prob, Storable prob) =>-      Int -> Trained typ sh prob -> ShowS--instance-   (Show typ, Shape.C sh, P.Show sh, P.Show prob, Storable prob) =>-      P.Show (T typ sh prob) where-   showsPrec = showsPrec--instance-   (Show typ, Shape.C sh, P.Show sh, P.Show prob, Storable prob) =>-      P.Show (Trained typ sh prob) where-   showsPrec = showsPrecTrained---class NFData typ where-   rnf ::-      (DeepSeq.NFData sh, DeepSeq.NFData prob, Shape.C sh) =>-      T typ sh prob -> ()-   rnfTrained ::-      (DeepSeq.NFData sh, DeepSeq.NFData prob, Shape.C sh) =>-      Trained typ sh prob -> ()--instance-   (NFData typ, DeepSeq.NFData sh, DeepSeq.NFData prob, Shape.C sh) =>-      DeepSeq.NFData (T typ sh prob) where-   rnf = rnf--instance-   (NFData typ, DeepSeq.NFData sh, DeepSeq.NFData prob, Shape.C sh) =>-      DeepSeq.NFData (Trained typ sh prob) where-   rnf = rnfTrained---class Format typ where-   format ::-      (Shape.C sh, Output out, Class.Real prob) =>-      String -> T typ sh prob -> out--instance-   (Format typ, Shape.C sh, Class.Real prob) =>-      Format.Format (T typ sh prob) where-   format = format----class Info typ where-   statesShape :: (Shape.C sh) => T typ sh prob -> sh-   statesShapeTrained :: (Shape.C sh) => Trained typ sh prob -> sh--class Generate typ where-   generate ::-      (Shape.Indexed sh, Class.Real prob, Rnd.Random prob, Rnd.RandomGen g) =>-      T typ sh prob -> Shape.Index sh -> MS.State g (Emission typ prob)--class EmissionProb typ where-   mapStatesShape ::-      (Shape.C sh0, Shape.C sh1) =>-      (sh0 -> sh1) -> T typ sh0 prob -> T typ sh1 prob-   {--   This function could be implemented generically in terms of emissionStateProb-   but that would require an Info constraint.-   -}-   emissionProb ::-      (Shape.C sh, Class.Real prob) =>-      T typ sh prob -> Emission typ prob -> Vector sh prob-   emissionStateProb ::-      (Shape.Indexed sh, Class.Real prob) =>-      T typ sh prob -> Emission typ prob -> Shape.Index sh -> prob-   emissionStateProb distr e s = emissionProb distr e StorableArray.! s--class (EmissionProb typ) => Estimate typ where-   accumulateEmissions ::-      (Shape.Indexed sh, Class.Real prob, Shape.Index sh ~ state) =>-      sh -> NonEmpty.T [] (state, Emission typ prob) -> Trained typ sh prob-   trainVector ::-      (Shape.C sh, Eq sh, Class.Real prob) =>-      Emission typ prob -> Vector sh prob -> Trained typ sh prob-   combine ::-      (Shape.C sh, Eq sh, Class.Real prob) =>-      Trained typ sh prob -> Trained typ sh prob -> Trained typ sh prob-   normalize ::-      (Shape.C sh, Eq sh, Class.Real prob) =>-      Trained typ sh prob -> T typ sh prob--accumulateEmissionVectors ::-   (Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>-   NonEmpty.T [] (Emission typ prob, Vector sh prob) -> Trained typ sh prob-accumulateEmissionVectors = NonEmpty.foldl1Map combine (uncurry trainVector)--instance-   (Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>-      Sg.Semigroup (Trained typ sh prob) where-   (<>) = combine---data Discrete symbol--newtype instance T (Discrete symbol) sh prob =-      Discrete (Matrix.General (Set symbol) sh prob)--newtype instance Trained (Discrete symbol) sh prob =-      DiscreteTrained (NonEmptyMap.T symbol (Vector sh prob))--type instance Emission (Discrete symbol) prob = symbol---instance (P.Show symbol, Ord symbol) => Show (Discrete symbol) where-   showsPrec prec (Discrete m) = P.showsPrec prec m-   showsPrecTrained prec (DiscreteTrained m) = P.showsPrec prec m--instance (DeepSeq.NFData symbol) => NFData (Discrete symbol) where-   rnf (Discrete m) = DeepSeq.rnf m-   rnfTrained (DiscreteTrained m) = DeepSeq.rnf m--instance (P.Show symbol, Ord symbol) => Format (Discrete symbol) where-   format fmt (Discrete m) =-      Output.formatAligned $-      map (\(sym,v) ->-            map (Identity . Output.text) $-            (show sym ++ ":") : map (printFmt fmt) (Vector.toList v)) $-      Array.toAssociations $ Matrix.toRowArray m---- cf. Data.Bifunctor.Flip-newtype Flip f b a = Flip {getFlip :: f a b}--printFmt :: (Class.Real a) => String -> a -> String-printFmt fmt =-   getFlip $ Class.switchReal (Flip $ printf fmt) (Flip $ printf fmt)--instance (Ord symbol) => Info (Discrete symbol) where-   statesShape (Discrete m) = Matrix.width m-   statesShapeTrained (DiscreteTrained m) = discreteStateShape m--instance (Ord symbol) => Generate (Discrete symbol) where-   generate (Discrete m) =-      randomItemProp . StorableArray.toAssociations . Matrix.takeColumn m--instance (Ord symbol) => EmissionProb (Discrete symbol) where-   mapStatesShape f (Discrete m) = Discrete $ Matrix.mapWidth f m-   emissionProb (Discrete m) = Matrix.takeRow m-   emissionStateProb (Discrete m) x s = m #! (x,s)--instance (Ord symbol) => Estimate (Discrete symbol) where-   accumulateEmissions sh =-      DiscreteTrained .-      NonEmptyMap.map-         (StorableArray.reshape sh .-          StorableArray.fromAssociations 0 (Shape.Deferred sh) .-          Map.toList) .-      NonEmptyMap.fromListWith (Map.unionWith (+)) .-      fmap (\(state,sym) -> (sym, Map.singleton (Shape.deferIndex sh state) 1))-   trainVector sym = DiscreteTrained . NonEmptyMap.singleton sym-   combine (DiscreteTrained distr0) (DiscreteTrained distr1) =-      DiscreteTrained $ NonEmptyMap.unionWith Vector.add distr0 distr1-   normalize (DiscreteTrained distr) =-      Discrete $ normalizeProbColumns $ discreteFromMap distr--normalizeProbColumns ::-   (Shape.C height, Shape.C width, Eq width, Class.Real a) =>-   Matrix.General height width a -> Matrix.General height width a-normalizeProbColumns m = m #/\ Matrix.columnSums m--discreteStateShape ::-   (Shape.C sh) => NonEmptyMap.T symbol (Vector sh prob) -> sh-discreteStateShape =-   StorableArray.shape . snd . fst . NonEmptyMap.minViewWithKey--discreteFromMap ::-   (Ord symbol, Shape.C sh, Eq sh, Class.Real prob) =>-   NonEmptyMap.T symbol (Vector sh prob) -> Matrix.General (Set symbol) sh prob-discreteFromMap m =-   Matrix.fromRowArray (discreteStateShape m) $-   Array.fromMap $ NonEmptyMap.flatten m--discreteFromList ::-   (Ord symbol, Shape.C sh, Eq sh, Class.Real prob) =>-   NonEmpty.T [] (symbol, Vector sh prob) -> T (Discrete symbol) sh prob-discreteFromList = Discrete . discreteFromMap . NonEmptyMap.fromList----data Gaussian emiSh--newtype instance T (Gaussian emiSh) stateSh a =-   Gaussian (Array stateSh (a, Vector emiSh a, Triangular.Upper emiSh a))--newtype instance Trained (Gaussian emiSh) stateSh a =-   GaussianTrained-      (StorableArray.Array (stateSh, MatrixShape.Hermitian (():+:emiSh)) a)--type instance Emission (Gaussian emiSh) a = Vector emiSh a---instance (Shape.C emiSh, P.Show emiSh) => Show (Gaussian emiSh) where-   showsPrec prec (Gaussian m) = P.showsPrec prec m-   showsPrecTrained prec (GaussianTrained m) = P.showsPrec prec m--instance (DeepSeq.NFData emiSh) => NFData (Gaussian emiSh) where-   rnf (Gaussian params) = DeepSeq.rnf params-   rnfTrained (GaussianTrained params) = DeepSeq.rnf params---instance (FormatArray emiSh) => Format (Gaussian emiSh) where-   format = runFormatGaussian $ Class.switchReal formatGaussian formatGaussian--newtype FormatGaussian out emiSh stateSh a =-   FormatGaussian-      {runFormatGaussian :: String -> T (Gaussian emiSh) stateSh a -> out}--formatGaussian ::-   (FormatArray emiSh, Shape.C stateSh,-    Class.Real a, Format.Format a, Output out) =>-   FormatGaussian out emiSh stateSh a-formatGaussian =-   FormatGaussian $ \fmt (Gaussian params) ->-      Format.format fmt $ Array.toList params---instance Info (Gaussian emiSh) where-   statesShape (Gaussian params) = Array.shape params-   statesShapeTrained (GaussianTrained params) =-      fst $ StorableArray.shape params--instance (Shape.C emiSh, Eq emiSh) => Generate (Gaussian emiSh) where-   generate (Gaussian allParams) state = do-      let (_c, center, covarianceChol) = 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) => EmissionProb (Gaussian emiSh) where-   mapStatesShape f (Gaussian m) = Gaussian $ Array.mapShape f m-   emissionProb (Gaussian allParams) x =-      StorableArray.fromBoxed $ fmap (gaussianEmissionProb x) allParams-   emissionStateProb (Gaussian allParams) x s =-      gaussianEmissionProb x $ allParams ! s--gaussianEmissionProb ::-   (Shape.C emiSh, Eq emiSh, Class.Real a) =>-   Vector emiSh a -> (a, Vector emiSh a, Triangular.Upper emiSh a) -> a-gaussianEmissionProb x (c, center, covarianceChol) =-   c * expSquared (Vector.sub x center -/# covarianceChol)--expSquared :: (Shape.C sh, Class.Real a) => Vector sh a -> a-expSquared =-   getNorm $ Class.switchReal (Norm expSquaredAux) (Norm expSquaredAux)--newtype Norm f a = Norm {getNorm :: f a -> a}--expSquaredAux ::-   (Shape.C sh, Class.Floating a, Vector.RealOf a ~ ar, Class.Real ar) =>-   Vector sh a -> ar-expSquaredAux x = exp ((-1/2) * Vector.norm2Squared x)---instance (Shape.C emiSh, Eq emiSh) => Estimate (Gaussian emiSh) where-   accumulateEmissions sh xs =-      let emiSh = StorableArray.shape $ snd $ NonEmpty.head xs-          hermSh = MatrixShape.hermitian MatrixShape.RowMajor (():+:emiSh)-      in GaussianTrained $-         Matrix.toRowMajor . Matrix.fromRowArray hermSh . Array.reshape sh .-         Array.accumulate Vector.add-            (Array.replicate (Shape.Deferred sh) (Vector.zero hermSh)) .-         map (\(state,v) -> (Shape.deferIndex sh state, extendedHermitian v)) .-         NonEmpty.flatten-            $ xs-   trainVector xs probs =-      GaussianTrained $ Matrix.toRowMajor $ probs |*- extendedHermitian xs-   combine (GaussianTrained m0) (GaussianTrained m1) =-      GaussianTrained $ Vector.add m0 m1-   {--     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 m) =-      let params (weight, centerSum, covarianceSum) =-             let c = recip (weight#!((),()))-                 center = Vector.scale c $ Matrix.flattenRow centerSum-             in  (center,-                  Matrix.sub-                     (Matrix.scaleRealReal c covarianceSum)-                     (Hermitian.outer MatrixShape.RowMajor center))-      in Gaussian $-         fmap (gaussianParameters . params .-               Hermitian.split . ArrMatrix.fromVector) $-         Matrix.toRowArray $ Matrix.fromRowMajor m--extendedHermitian ::-   (Shape.C emiSh, Class.Floating a) =>-   StorableArray.Array emiSh a ->-   StorableArray.Array (MatrixShape.Hermitian (():+:emiSh)) a-extendedHermitian =-   ArrMatrix.toVector .-   Hermitian.outer MatrixShape.RowMajor . Vector.append (Vector.one ())--{- |-input array must be non-empty--}-gaussianTrained ::-   (Shape.C emiSh, Eq emiSh, Shape.C stateSh, Class.Real prob) =>-   Array stateSh (prob, Vector emiSh prob, Matrix.Hermitian emiSh prob) ->-   Trained (Gaussian emiSh) stateSh prob-gaussianTrained =-   GaussianTrained . Matrix.toRowMajor .-   matrixFromRowArray "HMM.Distribution.gaussianTrained" .-   fmap-      (\(weight, center, covariance) ->-         ArrMatrix.toVector $-         Hermitian.stack-            (Hermitian.fromList MatrixShape.RowMajor () [weight])-            (Matrix.singleRow MatrixShape.RowMajor center)-            covariance)--matrixFromRowArray ::-   (Shape.C width, Eq width, Shape.C height, Class.Real a) =>-   String ->-   Array height (StorableArray.Array width a) ->-   Matrix.General height width a-matrixFromRowArray name xs =-   case Array.toList xs of-      [] -> error $ name ++ ": empty array"-      x:_ -> Matrix.fromRowArray (StorableArray.shape x) xs--gaussian ::-   (Shape.C emiSh, Shape.C stateSh, Class.Real prob) =>-   Array stateSh (Vector emiSh prob, Matrix.Hermitian emiSh prob) ->-   T (Gaussian emiSh) stateSh prob-gaussian = Gaussian . fmap gaussianParameters--gaussianParameters ::-   (Shape.C emiSh, Class.Real prob) =>-   (Vector emiSh prob, Matrix.Hermitian emiSh prob) ->-   (prob, Vector emiSh prob, Triangular.Upper emiSh prob)-gaussianParameters (center, covariance) =-   gaussianFromCholesky center $ HermitianPD.decompose covariance--gaussianFromCholesky ::-   (Shape.C emiSh, Class.Real prob) =>-   Vector emiSh prob -> Triangular.Upper emiSh prob ->-   (prob, Vector emiSh prob, Triangular.Upper emiSh prob)-gaussianFromCholesky center covarianceChol =-   let covarianceSqrtDet =-         Vector.product $ Triangular.takeDiagonal covarianceChol-   in  (recip (sqrt2pi ^ vectorDim center * covarianceSqrtDet),-        center, covarianceChol)--sqrt2pi :: (Class.Real a) => a-sqrt2pi = runIdentity $ Class.switchReal sqrt2piAux sqrt2piAux--sqrt2piAux :: (Floating a) => Identity a-sqrt2piAux = Identity $ sqrt (2*pi)---class ToCSV typ where-   toCells ::-      (Shape.C sh, Class.Real prob, P.Show prob) =>-      T typ sh prob -> [[String]]--class FromCSV typ where-   parseCells ::-      (Shape.C sh, Eq sh, Class.Real prob, Read prob) =>-      sh -> HMMCSV.CSVParser (T typ sh prob)--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 (CSVSymbol symbol) => ToCSV (Discrete symbol) where-   toCells (Discrete m) =-      map-         (\(symbol, probs) ->-            cellFromSymbol symbol : HMMCSV.cellsFromVector probs) $-      Array.toAssociations $ Matrix.toRowArray m--instance (CSVSymbol symbol) => FromCSV (Discrete symbol) where-   parseCells n =-      let p = parseSymbolProb n-      in fmap discreteFromList $-         liftA2 NonEmpty.Cons (HMMCSV.getRow >>= p) (HMMCSV.manyRowsUntilEnd p)--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) => ToCSV (Gaussian emiSh) where-   toCells (Gaussian params) =-      List.intercalate [[]] $-      map-         (\(_, center, covarianceChol) ->-            HMMCSV.cellsFromVector center :-            HMMCSV.cellsFromSquare (Triangular.toSquare covarianceChol)) $-      Array.toList params--instance (emiSh ~ Matrix.ShapeInt) => FromCSV (Gaussian emiSh) 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 . snd3) gs-      HMMCSV.assert (ListHT.allEqual sizes) $-         printf "dimensions of emissions mismatch: %s" (show sizes)-      return $ Gaussian $ Array.fromList sh gs--parseSingleGaussian ::-   (emiSh ~ Matrix.ShapeInt, Class.Real prob, Eq prob, Read prob) =>-   HMMCSV.CSVParser (prob, Vector emiSh prob, Triangular.Upper emiSh 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) =>-   Matrix.Square sh a -> Triangular.Upper sh a -> Bool-isUpperTriang m mt =-   Vector.toList (ArrMatrix.toVector m)-   ==-   Vector.toList (ArrMatrix.toVector (Triangular.toSquare mt))+import Math.HiddenMarkovModel.Public.Distribution

@@ -1,123 +0,0 @@-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.one stateSet,-      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

@@ -6,85 +6,7 @@ -} module Math.HiddenMarkovModel.Example.SineWave {-# WARNING "do not import that module, it is only intended for demonstration" #-}+   (module Math.HiddenMarkovModel.Example.SineWavePrivate)    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, singleton)--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.one stateSet,-      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+import Math.HiddenMarkovModel.Example.SineWavePrivate

@@ -1,157 +0,0 @@-{-# 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 qualified Data.Array.Comfort.Shape as Shape--import Text.Read.HT (maybeRead)--import Control.DeepSeq (NFData(rnf))-import Control.Monad (liftM2)--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 _ = ()--{- |-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 | 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.discreteFromList $-            (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.discreteFromList $-            (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)

@@ -1,157 +0,0 @@-{-# 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, revealGen, sumTransitions)-import Math.HiddenMarkovModel.Utility (normalizeFactor, normalizeProb)--import qualified Numeric.LAPACK.Matrix as Matrix-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.Shape as Shape--import qualified Data.NonEmpty.Class as NonEmptyC-import qualified Data.NonEmpty as NonEmpty-import qualified Data.Foldable as Fold-import Data.Traversable (Traversable)---{- |-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 typ, Shape.C sh, Eq sh, Floating prob,-    Class.Real prob, Distr.Emission typ prob ~ emission,-    Traversable f) =>-   T typ sh prob -> NonEmpty.T f emission -> prob-logLikelihood hmm = Fold.sum . fmap (log . fst) . alpha hmm--alpha ::-   (Distr.EmissionProb typ, Shape.C sh, Eq sh,-    Class.Real prob, Distr.Emission typ prob ~ emission,-    Traversable f) =>-   T typ 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 typ, Shape.C sh, Eq sh,-    Class.Real prob, Distr.Emission typ prob ~ emission,-    Traversable f, NonEmptyC.Reverse f) =>-   T typ 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.one $ StorableArray.shape $ initial hmm)--alphaBeta ::-   (Distr.EmissionProb typ, Shape.C sh, Eq sh,-    Class.Real prob, Distr.Emission typ prob ~ emission,-    Traversable f, NonEmptyC.Zip f, NonEmptyC.Reverse f) =>-   T typ 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 typ, Shape.C sh, Eq sh,-    Class.Real prob, Distr.Emission typ prob ~ emission,-    Traversable f, NonEmptyC.Zip f) =>-   T typ sh prob ->-   NonEmpty.T f emission ->-   NonEmpty.T f (prob, Vector sh prob) ->-   NonEmpty.T f (Vector sh prob) ->-   f (Matrix.Square sh prob)-xiFromAlphaBeta hmm xs calphas betas =-   let (cs,alphas) = Functor.unzip calphas-   in  NonEmptyC.zipWith4-          (\x alpha0 c1 beta1 ->-             Matrix.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 typ, Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ state,-    Distr.Emission typ prob ~ emission, Class.Real prob, Traversable f) =>-   T typ sh prob -> NonEmpty.T f emission -> NonEmpty.T f state-reveal = revealGen normalizeProb---{- |-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 typ, Shape.C sh, Eq sh,-    Class.Real prob, Distr.Emission typ prob ~ emission) =>-   T typ sh prob -> NonEmpty.T [] emission -> Trained typ 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.accumulateEmissionVectors $ NonEmptyC.zip xs zetas-       }

@@ -35,11 +35,11 @@ import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix import qualified Numeric.LAPACK.Matrix as Matrix 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.Static as ShapeStatic import qualified Data.Array.Comfort.Shape as Shape  import qualified Data.FixedLength as FL

@@ -1,331 +0,0 @@-{-# LANGUAGE TypeFamilies #-}-module Math.HiddenMarkovModel.Private where--import qualified Math.HiddenMarkovModel.Distribution as Distr-import qualified Math.HiddenMarkovModel.CSV as HMMCSV-import Math.HiddenMarkovModel.Utility (diagonal)--import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix-import qualified Numeric.LAPACK.Matrix.Square as Square-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.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.Semigroup ((<>))-import Data.Traversable (Traversable, mapAccumL)-import Data.Tuple.HT (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\".--}-data T typ sh prob =-   Cons {-      initial :: Vector sh prob,-      transition :: Matrix.Square sh prob,-      distribution :: Distr.T typ sh prob-   }-   deriving (Show)--instance-   (Distr.NFData typ, NFData sh, Shape.C sh, NFData prob, Storable prob) =>-      NFData (T typ sh prob) where-   rnf (Cons initial_ transition_ distribution_) =-      rnf (initial_, transition_, distribution_)--instance-   (Distr.Format typ, Format.FormatArray sh, Class.Real prob) =>-      Format.Format (T typ sh prob) where-   format fmt (Cons initial_ transition_ distribution_) =-      Format.format fmt (initial_, transition_, distribution_)--mapStatesShape ::-   (Distr.EmissionProb typ, Shape.C sh0, Shape.C sh1) =>-   (sh0 -> sh1) -> T typ sh0 prob -> T typ sh1 prob-mapStatesShape f hmm =-   Cons {-      initial = StorableArray.mapShape f $ initial hmm,-      transition = Square.mapSize f $ transition hmm,-      distribution = Distr.mapStatesShape f $ distribution hmm-   }---emission ::-   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob) =>-   T typ sh prob -> Distr.Emission typ prob -> Vector sh prob-emission  =  Distr.emissionProb . distribution---forward ::-   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,-    Distr.Emission typ prob ~ emission, Traversable f) =>-   T typ sh prob -> NonEmpty.T f emission -> prob-forward hmm = Vector.sum . NonEmpty.last . alpha hmm--alpha ::-   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,-    Distr.Emission typ prob ~ emission, Traversable f) =>-   T typ 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 ::-   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,-    Distr.Emission typ prob ~ emission, Traversable f) =>-   T typ sh prob -> NonEmpty.T f emission -> prob-backward hmm (NonEmpty.Cons x xs) =-   Vector.dot (initial hmm) $-   Vector.mul (emission hmm x) $-   NonEmpty.head $ beta hmm xs--beta ::-   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,-    Distr.Emission typ prob ~ emission, Traversable f) =>-   T typ 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.one $ StorableArray.shape $ initial hmm)---alphaBeta ::-   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,-    Distr.Emission typ prob ~ emission, Traversable f) =>-   T typ 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 ::-   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob) =>-   T typ sh prob -> Distr.Emission typ prob ->-   Vector sh prob -> Vector sh prob -> Matrix.Square sh prob-biscaleTransition hmm x alpha0 beta1 =-   (diagonal (Vector.mul (emission hmm x) beta1)-    #*##-    transition hmm)-   ##*#-   diagonal alpha0--xiFromAlphaBeta ::-   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,-    Distr.Emission typ prob ~ emission) =>-   T typ sh prob -> prob ->-   NonEmpty.T [] emission ->-   NonEmpty.T [] (Vector sh prob) ->-   NonEmpty.T [] (Vector sh prob) ->-   [Matrix.Square sh prob]-xiFromAlphaBeta hmm recipLikelihood xs alphas betas =-   zipWith3-      (\x alpha0 beta1 ->-         Matrix.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) =>-   [Matrix.Square 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 ::-   (Distr.EmissionProb typ, Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ state,-    Distr.Emission typ prob ~ emission, Class.Real prob, Traversable f) =>-   T typ sh prob -> NonEmpty.T f emission -> NonEmpty.T f state-reveal = revealGen id--revealGen ::-   (Distr.EmissionProb typ, Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ state,-    Distr.Emission typ prob ~ emission, Class.Real prob, Traversable f) =>-   (Vector (Shape.Deferred sh) prob -> Vector (Shape.Deferred sh) prob) ->-   T typ sh prob -> NonEmpty.T f emission -> NonEmpty.T f state-revealGen normalize hmm =-   fmap (Shape.revealIndex (StorableArray.shape $ initial hmm)) .-   revealStorable normalize (mapStatesShape Shape.Deferred hmm)--revealStorable ::-   (Distr.EmissionProb typ, Shape.InvIndexed sh, Eq sh,-    Shape.Index sh ~ state, Storable state,-    Distr.Emission typ prob ~ emission, Class.Real prob, Traversable f) =>-   (Vector sh prob -> Vector sh prob) ->-   T typ sh prob -> NonEmpty.T f emission -> NonEmpty.T f state-revealStorable normalize hmm (NonEmpty.Cons x xs) =-   uncurry (NonEmpty.scanr (StorableArray.!)) $-   mapFst (fst . Vector.argAbsMaximum) $-   mapAccumL-      (\alphai xi ->-         swap $ mapSnd (Vector.mul (emission hmm xi)) $-         matrixMaxMul (transition hmm) $ normalize alphai)-      (Vector.mul (emission hmm x) (initial hmm)) xs--matrixMaxMul ::-   (Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ ix, Storable ix,-    Class.Real a) =>-   Matrix.Square sh a -> Vector sh a ->-   (Vector sh ix, Vector sh a)-matrixMaxMul m v = Matrix.rowArgAbsMaximums $ Matrix.scaleColumns v 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 typ sh prob =-   Trained {-      trainedInitial :: Vector sh prob,-      trainedTransition :: Matrix.Square sh prob,-      trainedDistribution :: Distr.Trained typ sh prob-   }-   deriving (Show)--instance-   (Distr.NFData typ, NFData sh, Shape.C sh, NFData prob, Storable prob) =>-      NFData (Trained typ sh prob) where-   rnf hmm =-      rnf (trainedInitial hmm, trainedTransition hmm, trainedDistribution hmm)---sumTransitions ::-   (Shape.C sh, Eq sh, Class.Real e) =>-   T typ sh e -> [Matrix.Square sh e] -> Matrix.Square sh e-sumTransitions hmm =-   List.foldl' Matrix.add $-   Matrix.zero $ ArrMatrix.shape $ transition hmm--{- |-Baum-Welch algorithm--}-trainUnsupervised ::-   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob,-    Distr.Emission typ prob ~ emission) =>-   T typ sh prob -> NonEmpty.T [] emission -> Trained typ 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.accumulateEmissionVectors $ NonEmptyC.zip xs zetas-       }---mergeTrained ::-   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>-   Trained typ sh prob -> Trained typ sh prob -> Trained typ sh prob-mergeTrained hmm0 hmm1 =-   Trained {-      trainedInitial = Vector.add (trainedInitial hmm0) (trainedInitial hmm1),-      trainedTransition =-         Matrix.add (trainedTransition hmm0) (trainedTransition hmm1),-      trainedDistribution =-         trainedDistribution hmm0 <> trainedDistribution hmm1-   }--instance-   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>-      Sg.Semigroup (Trained typ sh prob) where-   (<>) = mergeTrained---toCells ::-   (Distr.ToCSV typ, Shape.Indexed sh, Class.Real prob, Show prob) =>-   T typ sh prob -> [[String]]-toCells hmm =-   (HMMCSV.cellsFromVector $ initial hmm) :-   (HMMCSV.cellsFromSquare $ transition hmm) ++-   [] :-   (Distr.toCells $ distribution hmm)--parseCSV ::-   (Distr.FromCSV typ, Shape.C stateSh, Eq stateSh,-    Class.Real prob, Read prob) =>-   (Int -> stateSh) -> HMMCSV.CSVParser (T typ 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-   }

@@ -1,259 +0,0 @@-{- |-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-         (squareFromLists, distance, matrixDistance)--import qualified Numeric.LAPACK.Matrix as Matrix-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 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.Map as NonEmptyMap-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.NonEmpty ((!:))-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.discreteFromList $-            ('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 =-   FL.curry-      (ShapeStatic.vector :: FL.T TypeNum.U4 Double -> Vector StateSet Double)---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 $-   case HMM.distribution hmm of Distr.Discrete m -> Matrix.takeRow m c--{- |-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)---{- |-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 :: ([Matrix.Square StateSet Double], [Matrix.Square 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 matrixDistance 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) ->-         (matrixDistance-             (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) ->-                (NonEmptyMap.size m0 == NonEmptyMap.size m1,-                 Fold.maximum $-                 Map.intersectionWith distance-                    (NonEmptyMap.flatten m0) (NonEmptyMap.flatten 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)) :-      []) ++-   []

@@ -1,88 +0,0 @@-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.Array as ArrMatrix-import qualified Numeric.LAPACK.Matrix as Matrix-import qualified Numeric.LAPACK.Vector as Vector-import Numeric.LAPACK.Matrix.Triangular (Diagonal)-import Numeric.LAPACK.Matrix.Array (ArrayMatrix)-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---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, 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---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 -> Matrix.Square sh a-squareConstant =-   (ArrMatrix.fromVector .) .-   Vector.constant . MatrixShape.square MatrixShape.RowMajor--squareFromLists ::-   (Shape.C sh, Eq sh, Storable a) => sh -> [Vector sh a] -> Matrix.Square 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---newtype Distance f a = Distance {getDistance :: f a -> f a -> a}--distance ::-   (Shape.C sh, Eq sh, Class.Real a) =>-   Vector sh a -> Vector sh a -> a-distance =-   getDistance $-   Class.switchReal-      (Distance $ (Vector.normInf .) . Vector.sub)-      (Distance $ (Vector.normInf .) . Vector.sub)--matrixDistance ::-   (Shape.C sh, Eq sh, Class.Real a) =>-   ArrayMatrix sh a -> ArrayMatrix sh a -> a-matrixDistance a b = distance (ArrMatrix.toVector a) (ArrMatrix.toVector b)

@@ -1,10 +1,15 @@ module Main where -import Math.HiddenMarkovModel.Test (tests)+import Test (tests) -import qualified Test.QuickCheck as QC+import qualified Test.DocTest.Driver as DocTest+import qualified Test.Main as TestMain   main :: IO ()-main =-   mapM_ (\(name,prop) -> putStr (name ++ ": ") >> QC.quickCheck prop) tests+main = DocTest.run $ do+   mapM_+      (\(name,prop) ->+         DocTest.printPrefix (name ++ ": ") >> DocTest.property prop)+      tests+   TestMain.main

@@ -0,0 +1,221 @@+module Test (tests) 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 Math.HiddenMarkovModel.Utility+         (squareFromLists, distance, matrixDistance)++import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Vector (Vector)++import qualified Data.Array.Comfort.Shape.Static as ShapeStatic+import qualified Data.Array.Comfort.Shape as Shape++import qualified Data.FixedLength as FL++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 qualified Data.NonEmpty.Class as NonEmptyC+import qualified Data.NonEmpty.Map as NonEmptyMap+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.NonEmpty ((!:))+++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.discreteFromList $+            ('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 =+   FL.curry+      (ShapeStatic.vector :: FL.T TypeNum.U4 Double -> Vector StateSet Double)+++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 $+   case HMM.distribution hmm of Distr.Discrete m -> Matrix.takeRow m c++{- |+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)+++{- |+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 :: ([Matrix.Square StateSet Double], [Matrix.Square 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 matrixDistance 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) ->+         (matrixDistance+             (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) ->+                (NonEmptyMap.size m0 == NonEmptyMap.size m1,+                 Fold.maximum $+                 Map.intersectionWith distance+                    (NonEmptyMap.flatten m0) (NonEmptyMap.flatten m1)))+++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) :+   []

@@ -0,0 +1,16 @@+-- Do not edit! Automatically created with doctest-extract.+module Test.Main where++import qualified Test.Math.HiddenMarkovModel.Example.TrafficLightPrivate+import qualified Test.Math.HiddenMarkovModel.Example.SineWavePrivate+import qualified Test.Math.HiddenMarkovModel.Example.CirclePrivate+import qualified Test.Math.HiddenMarkovModel.Normalized++import qualified Test.DocTest.Driver as DocTest++main :: DocTest.T ()+main = do+    Test.Math.HiddenMarkovModel.Example.TrafficLightPrivate.test+    Test.Math.HiddenMarkovModel.Example.SineWavePrivate.test+    Test.Math.HiddenMarkovModel.Example.CirclePrivate.test+    Test.Math.HiddenMarkovModel.Normalized.test

@@ -0,0 +1,59 @@+-- Do not edit! Automatically created with doctest-extract from private/Math/HiddenMarkovModel/Example/CirclePrivate.hs+{-# LINE 27 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}++module Test.Math.HiddenMarkovModel.Example.CirclePrivate where++import Math.HiddenMarkovModel.Example.CirclePrivate+import Test.DocTest.Base+import qualified Test.DocTest.Driver as DocTest++{-# LINE 28 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+import     qualified Math.HiddenMarkovModel as HMM+import     qualified Data.NonEmpty as NonEmpty+import     Data.Eq.HT (equating)++checkTraining     :: (Int, HMM) -> Bool+checkTraining     (maxDiff,hmm_) =+       maxDiff >=+       (length $ filter id $ NonEmpty.flatten $+        NonEmpty.zipWith (/=)+          (HMM.reveal hmm_ circle) (fmap fst circleLabeled))++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.CirclePrivate:61: "+{-# LINE 61 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+ DocTest.property+{-# LINE 61 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+     (checkTraining (0, hmm))+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.CirclePrivate:107: "+{-# LINE 107 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+ DocTest.property+{-# LINE 107 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+     (equating (take 1000 . NonEmpty.flatten) revealed $ fmap fst circleLabeled)+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.CirclePrivate:104: "+{-# LINE 104 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+ DocTest.example+{-# LINE 104 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+   (take 20 $ NonEmpty.flatten revealed)+  [ExpectedLine [LineChunk "[Q1,Q1,Q1,Q1,Q2,Q2,Q2,Q3,Q3,Q3,Q4,Q4,Q4,Q1,Q1,Q1,Q2,Q2,Q2,Q3]"]]+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.CirclePrivate:127: "+{-# LINE 127 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+ DocTest.property+{-# LINE 127 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+     (checkTraining (0, reconstructModel))+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.CirclePrivate:140: "+{-# LINE 140 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+ DocTest.property+{-# LINE 140 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+     (checkTraining (0, hmmTrainedSupervised))+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.CirclePrivate:147: "+{-# LINE 147 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+ DocTest.property+{-# LINE 147 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+     (checkTraining (0, hmmTrainedUnsupervised))+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.CirclePrivate:154: "+{-# LINE 154 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+ DocTest.property+{-# LINE 154 "private/Math/HiddenMarkovModel/Example/CirclePrivate.hs" #-}+     (checkTraining (40, hmmIterativelyTrained))

@@ -0,0 +1,26 @@+-- Do not edit! Automatically created with doctest-extract from private/Math/HiddenMarkovModel/Example/SineWavePrivate.hs+{-# LINE 21 "private/Math/HiddenMarkovModel/Example/SineWavePrivate.hs" #-}++module Test.Math.HiddenMarkovModel.Example.SineWavePrivate where++import Math.HiddenMarkovModel.Example.SineWavePrivate+import Test.DocTest.Base+import qualified Test.DocTest.Driver as DocTest++{-# LINE 22 "private/Math/HiddenMarkovModel/Example/SineWavePrivate.hs" #-}+import     qualified Data.NonEmpty as NonEmpty++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.SineWavePrivate:61: "+{-# LINE 61 "private/Math/HiddenMarkovModel/Example/SineWavePrivate.hs" #-}+ DocTest.example+{-# LINE 61 "private/Math/HiddenMarkovModel/Example/SineWavePrivate.hs" #-}+   (take 20 $ map fst $ NonEmpty.flatten sineWaveLabeled)+  [ExpectedLine [LineChunk "[Rising,Rising,High,High,High,Falling,Falling,Falling,Low,Low,Low,Rising,Rising,Rising,Rising,High,High,High,Falling,Falling]"]]+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.SineWavePrivate:74: "+{-# LINE 74 "private/Math/HiddenMarkovModel/Example/SineWavePrivate.hs" #-}+ DocTest.example+{-# LINE 74 "private/Math/HiddenMarkovModel/Example/SineWavePrivate.hs" #-}+   (take 20 $ NonEmpty.flatten revealed)+  [ExpectedLine [LineChunk "[Rising,Rising,High,High,High,Falling,Falling,Falling,Low,Low,Low,Low,Rising,Rising,Rising,High,High,High,Falling,Falling]"]]

@@ -0,0 +1,48 @@+-- Do not edit! Automatically created with doctest-extract from private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs+{-# LINE 22 "private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs" #-}++module Test.Math.HiddenMarkovModel.Example.TrafficLightPrivate where++import Math.HiddenMarkovModel.Example.TrafficLightPrivate+import Test.DocTest.Base+import qualified Test.DocTest.Driver as DocTest++{-# LINE 23 "private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs" #-}+import     qualified Data.NonEmpty as NonEmpty+import     Control.DeepSeq (deepseq)++verifyRevelations     :: HMM -> [Bool]+verifyRevelations     hmm_ =+       map (verifyRevelation hmm_) (NonEmpty.flatten labeledSequences)++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.TrafficLightPrivate:61: "+{-# LINE 61 "private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs" #-}+ DocTest.example+{-# LINE 61 "private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs" #-}+   (verifyRevelations hmm)+  [ExpectedLine [LineChunk "[True,True]"]]+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.TrafficLightPrivate:85: "+{-# LINE 85 "private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs" #-}+ DocTest.example+{-# LINE 85 "private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs" #-}+   (verifyRevelations hmmDisturbed)+  [ExpectedLine [LineChunk "[True,True]"]]+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.TrafficLightPrivate:138: "+{-# LINE 138 "private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs" #-}+ DocTest.example+{-# LINE 138 "private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs" #-}+   (verifyRevelations hmmTrainedSupervised)+  [ExpectedLine [LineChunk "[True,True]"]]+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.TrafficLightPrivate:153: "+{-# LINE 153 "private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs" #-}+ DocTest.example+{-# LINE 153 "private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs" #-}+   (verifyRevelations hmmTrainedUnsupervised)+  [ExpectedLine [LineChunk "[True,True]"]]+ DocTest.printPrefix "Math.HiddenMarkovModel.Example.TrafficLightPrivate:163: "+{-# LINE 163 "private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs" #-}+ DocTest.property+{-# LINE 163 "private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs" #-}+     (deepseq hmmIterativelyTrained True)

@@ -0,0 +1,18 @@+-- Do not edit! Automatically created with doctest-extract from private/Math/HiddenMarkovModel/Normalized.hs+{-# LINE 34 "private/Math/HiddenMarkovModel/Normalized.hs" #-}++module Test.Math.HiddenMarkovModel.Normalized where++import Math.HiddenMarkovModel.Normalized+import qualified Test.DocTest.Driver as DocTest++{-# LINE 35 "private/Math/HiddenMarkovModel/Normalized.hs" #-}+import     qualified Data.NonEmpty as NonEmpty++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Math.HiddenMarkovModel.Normalized:133: "+{-# LINE 133 "private/Math/HiddenMarkovModel/Normalized.hs" #-}+ DocTest.property+{-# LINE 133 "private/Math/HiddenMarkovModel/Normalized.hs" #-}+     (\x xs -> nonEmptyScanr (-) x xs == NonEmpty.scanr (-) x (xs::[Int]))