diff --git a/Changes.md b/Changes.md new file mode 100644 --- /dev/null +++ b/Changes.md @@ -0,0 +1,3 @@ +## 0.1 + +* `Distribution.Estimate` turned into a multi-parameter type class. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2015, Henning Thielemann + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Henning Thielemann nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#! /usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/hmm-lapack.cabal b/hmm-lapack.cabal new file mode 100644 --- /dev/null +++ b/hmm-lapack.cabal @@ -0,0 +1,101 @@ +Name: hmm-lapack +Version: 0.3 +Synopsis: Hidden Markov Models using HMatrix primitives +Description: + Hidden Markov Models implemented using HMatrix data types and operations. + + . + It implements: + . + * generation of samples of emission sequences, + . + * computation of the likelihood of an observed sequence of emissions, + . + * construction of most likely state sequence + that produces an observed sequence of emissions, + . + * supervised and unsupervised training of the model by Baum-Welch algorithm. + . + It supports any kind of emission distribution, + where discrete and multivariate Gaussian distributions + are implemented as examples. + . + For an introduction please refer to the examples: + . + * "Math.HiddenMarkovModel.Example.TrafficLight" + . + * "Math.HiddenMarkovModel.Example.SineWave" + . + * "Math.HiddenMarkovModel.Example.Circle" + . + An alternative package without foreign calls is @hmm@. +Homepage: http://hub.darcs.net/thielema/hmm-hmatrix +License: BSD3 +License-File: LICENSE +Author: Henning Thielemann +Maintainer: haskell@henning-thielemann.de +Category: Math +Build-Type: Simple +Cabal-Version: >=1.10 +Extra-Source-Files: + Changes.md + +Source-Repository this + Tag: 0.3 + Type: darcs + Location: http://hub.darcs.net/thielema/hmm-hmatrix + +Source-Repository head + Type: darcs + Location: http://hub.darcs.net/thielema/hmm-hmatrix + +Library + Exposed-Modules: + Math.HiddenMarkovModel + Math.HiddenMarkovModel.Named + Math.HiddenMarkovModel.Distribution + Math.HiddenMarkovModel.Pattern + Math.HiddenMarkovModel.Example.TrafficLight + Math.HiddenMarkovModel.Example.SineWave + Math.HiddenMarkovModel.Example.Circle + Math.HiddenMarkovModel.Test + Other-Modules: + Math.HiddenMarkovModel.Example.TrafficLightPrivate + Math.HiddenMarkovModel.Example.CirclePrivate + Math.HiddenMarkovModel.Normalized + Math.HiddenMarkovModel.Private + Math.HiddenMarkovModel.Utility + Math.HiddenMarkovModel.CSV + Build-Depends: + lapack >=0.2 && <0.3, + fixed-length >=0.2 && <0.3, + tfp >=1.0 && <1.1, + netlib-ffi >=0.1.1 && <0.2, + comfort-array >=0.2 && <0.3, + QuickCheck >=2.5 && <3, + explicit-exception >=0.1.7 && <0.2, + boxes >=0.1.5 && <0.2, + lazy-csv >=0.5 && <0.6, + random >=1.0 && <1.2, + transformers >= 0.2 && <0.6, + non-empty >=0.2.1 && <0.4, + semigroups >=0.17 && <0.19, + containers >=0.4.2 && <0.7, + utility-ht >=0.0.12 && <0.1, + deepseq >=1.3 && <1.5, + prelude-compat >=0.0 && <0.1, + base >=4.5 && <5 + Hs-Source-Dirs: src + Default-Language: Haskell2010 + GHC-Options: -Wall + +Test-Suite hmm-test + Type: exitcode-stdio-1.0 + Build-Depends: + hmm-lapack, + QuickCheck, + base + Main-Is: Main.hs + Hs-Source-Dirs: test + Default-Language: Haskell2010 + GHC-Options: -Wall diff --git a/src/Math/HiddenMarkovModel.hs b/src/Math/HiddenMarkovModel.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel.hs @@ -0,0 +1,221 @@ +{-# LANGUAGE TypeFamilies #-} +module Math.HiddenMarkovModel ( + T(..), + Discrete, DiscreteTrained, + Gaussian, GaussianTrained, + uniform, + generate, + generateLabeled, + probabilitySequence, + Normalized.logLikelihood, + Normalized.reveal, + + Trained(..), + trainSupervised, + Normalized.trainUnsupervised, + mergeTrained, finishTraining, trainMany, + deviation, + + toCSV, + fromCSV, + ) where + +import qualified Math.HiddenMarkovModel.Distribution as Distr +import qualified Math.HiddenMarkovModel.Normalized as Normalized +import qualified Math.HiddenMarkovModel.CSV as HMMCSV +import Math.HiddenMarkovModel.Private + (T(..), Trained(..), mergeTrained, toCells, parseCSV) +import Math.HiddenMarkovModel.Utility + (SquareMatrix, squareConstant, + randomItemProp, normalizeProb, attachOnes) + +import qualified Numeric.LAPACK.Matrix as Matrix +import qualified Numeric.LAPACK.Vector as Vector +import qualified Numeric.LAPACK.Scalar as Scalar + +import qualified Numeric.Netlib.Class as Class + +import qualified Data.Array.Comfort.Storable as StorableArray +import qualified Data.Array.Comfort.Shape as Shape +import qualified Data.Array.Comfort.Boxed as Array + +import qualified Text.CSV.Lazy.String as CSV + +import qualified System.Random as Rnd + +import qualified Control.Monad.Exception.Synchronous as ME +import qualified Control.Monad.Trans.State as MS +import qualified Control.Monad.HT as Monad + +import qualified Data.NonEmpty as NonEmpty +import Data.Traversable (Traversable, mapAccumL) +import Data.Foldable (Foldable) + + + +type DiscreteTrained symbol sh prob = + Trained (Distr.DiscreteTrained symbol sh prob) sh prob +type Discrete symbol sh prob = T (Distr.Discrete symbol sh prob) sh prob + +type GaussianTrained emiSh stateSh a = + Trained (Distr.GaussianTrained emiSh stateSh a) stateSh a +type Gaussian emiSh stateSh a = T (Distr.Gaussian emiSh stateSh a) stateSh a + + +{- | +Create a model with uniform probabilities +for initial vector and transition matrix +given a distribution for the emissions. +You can use this as a starting point for 'Normalized.trainUnsupervised'. +-} +uniform :: + (Distr.Info distr, Distr.StateShape distr ~ sh, Shape.C sh, + Distr.Probability distr ~ prob) => + distr -> T distr sh prob +uniform distr = + let sh = Distr.statesShape distr + c = recip $ fromIntegral $ Shape.size sh + in Cons { + initial = Vector.constant sh c, + transition = squareConstant sh c, + distribution = distr + } + + +probabilitySequence :: + (Traversable f, Distr.EmissionProb distr, + Distr.StateShape distr ~ sh, Shape.Indexed sh, Shape.Index sh ~ state, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + T distr sh prob -> f (state, emission) -> f prob +probabilitySequence hmm = + snd + . + mapAccumL + (\index (s, e) -> + ((transition hmm StorableArray.!) . flip (,) s, + index s * Distr.emissionStateProb (distribution hmm) e s)) + (initial hmm StorableArray.!) + +generate :: + (Rnd.RandomGen g, Ord prob, Rnd.Random prob, Distr.Generate distr, + Distr.StateShape distr ~ sh, Shape.Indexed sh, Shape.Index sh ~ state, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + T distr sh prob -> g -> [emission] +generate hmm = map snd . generateLabeled hmm + +generateLabeled :: + (Rnd.RandomGen g, Ord prob, Rnd.Random prob, Distr.Generate distr, + Distr.StateShape distr ~ sh, Shape.Indexed sh, Shape.Index sh ~ state, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + T distr sh prob -> g -> [(state, emission)] +generateLabeled hmm = + MS.evalState $ + flip MS.evalStateT (initial hmm) $ + Monad.repeat $ MS.StateT $ \v0 -> do + s <- + randomItemProp $ + zip (Shape.indices $ StorableArray.shape v0) (Vector.toList v0) + x <- Distr.generate (distribution hmm) s + return ((s, x), Matrix.takeColumn (transition hmm) s) + + + +{- | +Contribute a manually labeled emission sequence to a HMM training. +-} +trainSupervised :: + (Distr.StateShape distr ~ sh, Shape.Index sh ~ state, + Distr.Estimate tdistr distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + sh -> NonEmpty.T [] (state, emission) -> Trained tdistr sh prob +trainSupervised sh xs = + let getState (s, _x) = s + in Trained { + trainedInitial = + StorableArray.fromAssociations sh 0 + [(getState (NonEmpty.head xs), 1)], + trainedTransition = + Matrix.transpose $ + StorableArray.accumulate (+) (squareConstant sh 0) $ + attachOnes $ NonEmpty.mapAdjacent (,) $ fmap getState xs, + trainedDistribution = + Distr.accumulateEmissions $ Array.map attachOnes $ + Array.accumulate (flip (:)) + (Array.fromList sh $ replicate (Shape.size sh) []) + (NonEmpty.flatten xs) + } + +finishTraining :: + (Shape.C sh, Eq sh, + Distr.Estimate tdistr distr, Distr.Probability distr ~ prob) => + Trained tdistr sh prob -> T distr sh prob +finishTraining hmm = + Cons { + initial = normalizeProb $ trainedInitial hmm, + transition = normalizeProbColumns $ trainedTransition hmm, + distribution = Distr.normalize $ trainedDistribution hmm + } + +normalizeProbColumns :: + (Shape.C sh, Eq sh, Class.Real a) => SquareMatrix sh a -> SquareMatrix sh a +normalizeProbColumns m = + Matrix.scaleColumns (StorableArray.map recip (Matrix.columnSums m)) m + +trainMany :: + (Shape.C sh, Eq sh, + Distr.Estimate tdistr distr, Distr.Probability distr ~ prob, + Foldable f) => + (trainingData -> Trained tdistr sh prob) -> + NonEmpty.T f trainingData -> T distr sh prob +trainMany train = + finishTraining . NonEmpty.foldl1Map mergeTrained train + + + + + +{- | +Compute maximum deviation between initial and transition probabilities. +You can use this as abort criterion for unsupervised training. +We omit computation of differences between the emission probabilities. +This simplifies matters a lot and +should suffice for defining an abort criterion. +-} +deviation :: + (Shape.InvIndexed sh, Eq sh, Class.Real prob, Ord prob) => + T distr sh prob -> T distr sh prob -> prob +deviation hmm0 hmm1 = + deviationVec (initial hmm0) (initial hmm1) + `max` + deviationVec (transition hmm0) (transition hmm1) + +deviationVec :: + (Shape.InvIndexed sh, Eq sh, Class.Real a) => + StorableArray.Array sh a -> StorableArray.Array sh a -> a +deviationVec = + getDeviation $ Class.switchReal deviationVecAux deviationVecAux + +newtype Deviation f a = Deviation {getDeviation :: f a -> f a -> a} + +deviationVecAux :: + (Shape.InvIndexed sh, Eq sh, Ord a, Class.Real a, Scalar.RealOf a ~ a) => + Deviation (StorableArray.Array sh) a +deviationVecAux = + Deviation $ \x y -> + Scalar.absolute $ snd $ Vector.argAbsMaximum $ Vector.sub x y + + +toCSV :: + (Distr.ToCSV distr, Shape.Indexed sh, Class.Real prob, Show prob) => + T distr sh prob -> String +toCSV hmm = + CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $ + toCells hmm + +fromCSV :: + (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh, + Shape.Indexed stateSh, Shape.Index stateSh ~ state, + Class.Real prob, Read prob) => + (Int -> stateSh) -> String -> ME.Exceptional String (T distr stateSh prob) +fromCSV makeShape = + MS.evalStateT (parseCSV makeShape) . map HMMCSV.fixShortRow . CSV.parseCSV diff --git a/src/Math/HiddenMarkovModel/CSV.hs b/src/Math/HiddenMarkovModel/CSV.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/CSV.hs @@ -0,0 +1,160 @@ +module Math.HiddenMarkovModel.CSV where + +import Math.HiddenMarkovModel.Utility (SquareMatrix, vectorDim) + +import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape +import qualified Numeric.LAPACK.Matrix as Matrix +import qualified Numeric.LAPACK.Vector as Vector +import Numeric.LAPACK.Matrix (ZeroInt) +import Numeric.LAPACK.Vector (Vector) + +import qualified Numeric.Netlib.Class as Class + +import qualified Data.Array.Comfort.Storable as ComfortArray +import qualified Data.Array.Comfort.Shape as Shape + +import qualified Text.CSV.Lazy.String as CSV +import Text.Read.HT (maybeRead) +import Text.Printf (printf) + +import qualified Control.Monad.Exception.Synchronous as ME +import qualified Control.Monad.Trans.Class as MT +import qualified Control.Monad.Trans.State as MS +import Control.Monad.Exception.Synchronous (Exceptional) +import Control.Monad (liftM2, replicateM, unless) + +import qualified Data.List.Reverse.StrictElement as Rev +import qualified Data.List.HT as ListHT + + +cellsFromVector :: + (Shape.C sh, Show a, Class.Real a) => Vector sh a -> [String] +cellsFromVector = map show . Vector.toList + +cellsFromSquare :: + (Shape.Indexed sh, Show a, Class.Real a) => SquareMatrix sh a -> [[String]] +cellsFromSquare = map (map show . Vector.toList) . Matrix.toRows + +padTable :: a -> [[a]] -> [[a]] +padTable x xs = + let width = maximum (map length xs) + in map (ListHT.padRight x width) xs + + +type CSVParser = MS.StateT CSV.CSVResult (Exceptional String) + +assert :: Bool -> String -> CSVParser () +assert cond msg = + unless cond $ MT.lift $ ME.throw msg + +retrieveShortRow :: CSV.CSVError -> Maybe CSV.CSVRow +retrieveShortRow err = + case err of + CSV.IncorrectRow {CSV.csvFields = row} -> Just row + _ -> Nothing + +fixShortRow :: + Either [CSV.CSVError] CSV.CSVRow -> Either [CSV.CSVError] CSV.CSVRow +fixShortRow erow = + case erow of + Left errs -> + case ListHT.partitionMaybe retrieveShortRow errs of + ([row], []) -> Right row + _ -> Left errs + _ -> erow + +maybeGetRow :: CSVParser (Maybe CSV.CSVRow) +maybeGetRow = do + csv0 <- MS.get + case csv0 of + [] -> return Nothing + item : csv1 -> do + MS.put csv1 + case item of + Right row -> return (Just row) + Left errors -> + MT.lift $ ME.throw $ unlines $ map CSV.ppCSVError errors + +getRow :: CSVParser CSV.CSVRow +getRow = + MT.lift . ME.fromMaybe "unexpected end of file" =<< maybeGetRow + +checkEmptyRow :: CSV.CSVRow -> Exceptional String () +checkEmptyRow row = + case filter (not . null . CSV.csvFieldContent) row of + [] -> return () + cell:_ -> ME.throw $ printf "%d: expected empty row" (CSV.csvRowNum cell) + +skipEmptyRow :: CSVParser () +skipEmptyRow = MT.lift . checkEmptyRow =<< getRow + +manySepUntilEnd :: CSVParser a -> CSVParser [a] +manySepUntilEnd p = + let go = liftM2 (:) p $ do + mrow <- maybeGetRow + case mrow of + Nothing -> return [] + Just row -> do + MT.lift $ checkEmptyRow row + go + in go + +manyRowsUntilEnd :: (CSV.CSVRow -> CSVParser a) -> CSVParser [a] +manyRowsUntilEnd p = + let go = do + mrow <- maybeGetRow + case mrow of + Nothing -> return [] + Just row -> liftM2 (:) (p row) go + in go + +parseVectorCells :: + (Read a, Class.Real a) => + CSVParser (Vector ZeroInt a) +parseVectorCells = + parseVectorFields =<< getRow + +-- ToDo: Maybe check row consistency already here? +parseVectorFields :: + (Read a, Class.Real a) => + CSV.CSVRow -> CSVParser (Vector ZeroInt a) +parseVectorFields = + MT.lift . fmap Vector.autoFromList . mapM parseNumberCell . + Rev.dropWhile (null . CSV.csvFieldContent) + +parseNonEmptyVectorCells :: + (Read a, Class.Real a) => + CSVParser (Vector ZeroInt a) +parseNonEmptyVectorCells = do + v <- parseVectorCells + assert (vectorDim v > 0) "no data for vector" + return v + +cellContent :: CSV.CSVField -> Exceptional String String +cellContent field = + case field of + CSV.CSVFieldError {} -> ME.throw $ CSV.ppCSVField field + CSV.CSVField { CSV.csvFieldContent = str } -> return str + +parseNumberCell :: (Read a) => CSV.CSVField -> Exceptional String a +parseNumberCell field = do + str <- cellContent field + ME.fromMaybe (printf "field content \"%s\" is not a number" str) $ + maybeRead str + +parseSquareMatrixCells :: + (Shape.C sh, Read a, Class.Real a) => + sh -> CSVParser (SquareMatrix sh a) +parseSquareMatrixCells sh = do + let n = Shape.size sh + rows <- replicateM n parseVectorCells + assert (not $ null rows) "no rows" + assert (all ((n==) . vectorDim) rows) "inconsistent matrix dimensions" + return $ + ComfortArray.reshape (MatrixShape.square MatrixShape.RowMajor sh) $ + Matrix.fromRows (Shape.ZeroBased n) rows + +parseStringList :: CSV.CSVRow -> CSVParser [String] +parseStringList = + MT.lift . mapM cellContent . + Rev.dropWhile (null . CSV.csvFieldContent) diff --git a/src/Math/HiddenMarkovModel/Distribution.hs b/src/Math/HiddenMarkovModel/Distribution.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Distribution.hs @@ -0,0 +1,481 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +module Math.HiddenMarkovModel.Distribution ( + Emission, Probability, StateShape, + Info(..), Generate(..), EmissionProb(..), Estimate(..), + + Discrete(..), DiscreteTrained(..), + Gaussian(..), GaussianTrained(..), gaussian, + + ToCSV(..), FromCSV(..), HMMCSV.CSVParser, CSVSymbol(..), + ) where + +import qualified Math.HiddenMarkovModel.CSV as HMMCSV +import Math.HiddenMarkovModel.Utility (SquareMatrix, randomItemProp, vectorDim) + +import qualified Numeric.LAPACK.Matrix.HermitianPositiveDefinite as HermitianPD +import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian +import qualified Numeric.LAPACK.Matrix.Triangular as Triangular +import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape +import qualified Numeric.LAPACK.Matrix as Matrix +import qualified Numeric.LAPACK.Vector as Vector +import Numeric.LAPACK.Matrix ((<#)) +import Numeric.LAPACK.Vector (Vector) +import Numeric.LAPACK.Format (FormatArray, Format(format)) + +import qualified Numeric.Netlib.Class as Class +import Foreign.Storable (Storable) + +import qualified Data.Array.Comfort.Storable as StorableArray +import qualified Data.Array.Comfort.Shape as Shape +import qualified Data.Array.Comfort.Boxed as Array +import Data.Array.Comfort.Boxed (Array, (!)) + +import qualified System.Random as Rnd + +import qualified Text.CSV.Lazy.String as CSV +import qualified Text.PrettyPrint.Boxes as TextBox +import Text.PrettyPrint.Boxes ((<>), (<+>)) +import Text.Read.HT (maybeRead) +import Text.Printf (printf) + +import qualified Control.Monad.Exception.Synchronous as ME +import qualified Control.Monad.Trans.Class as MT +import qualified Control.Monad.Trans.State as MS +import Control.DeepSeq (NFData, rnf) +import Control.Monad (liftM2) +import Control.Applicative (liftA2, (<|>)) + +import qualified Data.NonEmpty as NonEmpty +import qualified Data.Foldable as Fold +import qualified Data.Map as Map +import qualified Data.Set as Set +import qualified Data.List.HT as ListHT +import qualified Data.List as List +import Data.Functor.Identity (Identity(Identity), runIdentity) +import Data.Foldable (Foldable, foldMap) +import Data.Tuple.HT (mapFst, fst3, swap) +import Data.Monoid (Endo(Endo, appEndo)) +import Data.Map (Map) +import Data.Maybe (fromMaybe, listToMaybe) + +import Prelude () +import Prelude2010 + + +type HermitianMatrix sh = Hermitian.Hermitian sh +type UpperTriangular sh = Triangular.Upper sh + + +type family Probability distr +type family Emission distr +type family StateShape distr + + +class (Class.Real (Probability distr)) => Info distr where + statesShape :: distr -> StateShape distr + +class (Class.Real (Probability distr)) => Generate distr where + generate :: + (Rnd.RandomGen g, Emission distr ~ emission, StateShape distr ~ sh) => + distr -> Shape.Index sh -> MS.State g emission + +class + (Shape.Indexed (StateShape distr), Class.Real (Probability distr)) => + EmissionProb distr where + {- + This function could be implemented generically in terms of emissionStateProb + but that would require an Info constraint. + -} + emissionProb :: + distr -> Emission distr -> Vector (StateShape distr) (Probability distr) + emissionStateProb :: + distr -> Emission distr -> Shape.Index (StateShape distr) -> Probability distr + emissionStateProb distr e s = emissionProb distr e StorableArray.! s + +class + (Distribution tdistr ~ distr, Trained distr ~ tdistr, EmissionProb distr) => + Estimate tdistr distr where + type Distribution tdistr + type Trained distr + accumulateEmissions :: + (Probability distr ~ prob, StateShape distr ~ sh) => + Array sh [(Emission distr, prob)] -> tdistr + -- could as well be in Semigroup class + combine :: tdistr -> tdistr -> tdistr + normalize :: tdistr -> distr + + +newtype Discrete symbol sh prob = Discrete (Map symbol (Vector sh prob)) + deriving (Show) + +newtype + DiscreteTrained symbol sh prob = + DiscreteTrained (Map symbol (Vector sh prob)) + deriving (Show) + +type instance Probability (Discrete symbol sh prob) = prob +type instance Emission (Discrete symbol sh prob) = symbol +type instance StateShape (Discrete symbol sh prob) = sh + + +instance + (NFData sh, NFData prob, NFData symbol) => + NFData (Discrete symbol sh prob) where + rnf (Discrete m) = rnf m + +instance + (NFData sh, NFData prob, NFData symbol) => + NFData (DiscreteTrained symbol sh prob) where + rnf (DiscreteTrained m) = rnf m + +instance + (FormatArray sh, Class.Real prob, Format symbol) => + Format (Discrete symbol sh prob) where + format fmt (Discrete m) = + TextBox.vsep 1 TextBox.left $ + map (\(sym,v) -> format fmt sym <> TextBox.char ':' <+> format fmt v) $ + Map.toAscList m + +instance + (Shape.C sh, Class.Real prob, Ord symbol) => + Info (Discrete symbol sh prob) where + statesShape (Discrete m) = StorableArray.shape $ snd $ Map.findMin m + +instance + (Shape.Indexed sh, Class.Real prob, Ord symbol, Ord prob, Rnd.Random prob) => + Generate (Discrete symbol sh prob) where + generate (Discrete m) state = + randomItemProp $ Map.toAscList $ fmap (StorableArray.! state) m + +instance + (Shape.Indexed sh, Class.Real prob, Ord symbol) => + EmissionProb (Discrete symbol sh prob) where + emissionProb (Discrete m) = + mapLookup "emitDiscrete: unknown emission symbol" m + +instance + (Shape.Indexed sh, Eq sh, Class.Real prob, Ord symbol) => + Estimate (DiscreteTrained symbol sh prob) (Discrete symbol sh prob) where + type Distribution (DiscreteTrained symbol sh prob) = Discrete symbol sh prob + type Trained (Discrete symbol sh prob) = DiscreteTrained symbol sh prob + accumulateEmissions grouped = + let set = Set.toAscList $ foldMap (Set.fromList . map fst) grouped + emi = Map.fromAscList $ zip set [0..] + in DiscreteTrained $ Map.fromAscList $ zip set $ + transposeVectorList $ + Array.map + (StorableArray.accumulate (+) + (Vector.constant (Shape.ZeroBased $ length set) 0) . + map (mapFst + (mapLookup "estimateDiscrete: unknown emission symbol" emi))) + grouped + combine (DiscreteTrained distr0) (DiscreteTrained distr1) = + DiscreteTrained $ Map.unionWith Vector.add distr0 distr1 + normalize (DiscreteTrained distr) = + Discrete $ if Map.null distr then distr else normalizeProbVecs distr + +transposeVectorList :: + (Shape.C sh, Eq sh, Class.Real a) => + Array sh (Vector Matrix.ZeroInt a) -> [Vector sh a] +transposeVectorList xs = + case Array.toList xs of + [] -> [] + x:_ -> Matrix.toRows $ Matrix.fromColumnArray (StorableArray.shape x) xs + +normalizeProbVecs :: + (Shape.C sh, Eq sh, Foldable f, Functor f, Class.Real a) => + f (Vector sh a) -> f (Vector sh a) +normalizeProbVecs vs = + let factors = + StorableArray.map recip $ List.foldl1' Vector.add $ Fold.toList vs + in fmap (Vector.mul factors) vs + +mapLookup :: (Ord k) => String -> Map.Map k a -> k -> a +mapLookup msg dict x = Map.findWithDefault (error msg) x dict + + +newtype Gaussian emiSh stateSh a = + Gaussian (Array stateSh (Vector emiSh a, UpperTriangular emiSh a, a)) + deriving (Show) + +newtype GaussianTrained emiSh stateSh a = + GaussianTrained + (Array stateSh + (Maybe (Vector emiSh a, HermitianMatrix emiSh a, a))) + deriving (Show) + +type instance Probability (Gaussian emiSh stateSh a) = a +type instance Emission (Gaussian emiSh stateSh a) = Vector emiSh a +type instance StateShape (Gaussian emiSh stateSh a) = stateSh + + +instance + (NFData emiSh, NFData stateSh, Shape.C stateSh, NFData a, Storable a) => + NFData (Gaussian emiSh stateSh a) where + rnf (Gaussian params) = rnf params + +instance + (NFData emiSh, NFData stateSh, Shape.C stateSh, NFData a, Storable a) => + NFData (GaussianTrained emiSh stateSh a) where + rnf (GaussianTrained params) = rnf params + + +instance + (FormatArray emiSh, Shape.C stateSh, Class.Real a) => + Format (Gaussian emiSh stateSh a) where + format = runFormatGaussian $ Class.switchReal formatGaussian formatGaussian + +newtype FormatGaussian emiSh stateSh a = + FormatGaussian + {runFormatGaussian :: String -> Gaussian emiSh stateSh a -> TextBox.Box} + +formatGaussian :: + (FormatArray emiSh, Shape.C stateSh, Class.Real a, Format a) => + FormatGaussian emiSh stateSh a +formatGaussian = + FormatGaussian $ \fmt (Gaussian params) -> format fmt $ Array.toList params + + +instance + (Shape.Indexed stateSh, Eq stateSh, Class.Real a) => + Info (Gaussian emiSh stateSh a) where + statesShape (Gaussian params) = Array.shape params + +instance + (Shape.C emiSh, Eq emiSh, Shape.Indexed stateSh, Eq stateSh, Class.Real a) => + Generate (Gaussian emiSh stateSh a) where + generate (Gaussian allParams) state = do + let (center, covarianceChol, _c) = allParams ! state + seed <- MS.state Rnd.random + return $ + Vector.add center $ + Vector.random Vector.Normal (StorableArray.shape center) seed + <# covarianceChol + +instance + (Shape.C emiSh, Eq emiSh, Shape.Indexed stateSh, Eq stateSh, Class.Real a) => + EmissionProb (Gaussian emiSh stateSh a) where + emissionProb (Gaussian allParams) x = + Vector.fromList (Array.shape allParams) $ + map (emissionProbGen x) $ Array.toList allParams + emissionStateProb (Gaussian allParams) x s = + emissionProbGen x $ allParams ! s + +emissionProbGen :: + (Shape.C emiSh, Eq emiSh, Class.Real a) => + Vector emiSh a -> (Vector emiSh a, UpperTriangular emiSh a, a) -> a +emissionProbGen x (center, covarianceChol, c) = + let x0 = + Matrix.solveVector (Triangular.transpose covarianceChol) $ + Vector.sub x center + in c * cexp ((-1/2) * Vector.inner x0 x0) + + +instance + (Shape.C emiSh, Eq emiSh, Shape.Indexed stateSh, Eq stateSh, Class.Real a) => + Estimate + (GaussianTrained emiSh stateSh a) + (Gaussian emiSh stateSh a) where + type Distribution (GaussianTrained emiSh stateSh a) = + Gaussian emiSh stateSh a + type Trained (Gaussian emiSh stateSh a) = GaussianTrained emiSh stateSh a + accumulateEmissions = + let params xs = + (NonEmpty.foldl1Map Vector.add (uncurry $ flip Vector.scale) xs, + covarianceReal $ fmap swap xs, + Fold.sum $ fmap snd xs) + in GaussianTrained . fmap (fmap params . NonEmpty.fetch) + combine (GaussianTrained distr0) (GaussianTrained distr1) = + let comb (center0, covariance0, weight0) + (center1, covariance1, weight1) = + (Vector.add center0 center1, + Vector.add covariance0 covariance1, + weight0 + weight1) + in GaussianTrained $ Array.zipWith (maybePlus comb) distr0 distr1 + {- + Sum_i (xi-m) * (xi-m)^T + = Sum_i xi*xi^T + Sum_i m*m^T - Sum_i xi*m^T - Sum_i m*xi^T + = Sum_i xi*xi^T - Sum_i m*m^T + = Sum_i xi*xi^T - n * m*m^T + -} + normalize (GaussianTrained distr) = + let params (centerSum, covarianceSum, weight) = + let c = recip weight + center = Vector.scale c centerSum + in (center, + Vector.sub (Vector.scale c covarianceSum) + (Hermitian.outer MatrixShape.RowMajor center)) + in Gaussian $ + fmap + (gaussianParameters . params . + fromMaybe + (error "Distribution.normalize: undefined array element")) $ + distr + +-- ToDo: could be managed by semigroup +maybePlus :: (a -> a -> a) -> Maybe a -> Maybe a -> Maybe a +maybePlus f mx my = liftA2 f mx my <|> mx <|> my + + +newtype CovarianceReal f emiSh a = + CovarianceReal + {getCovarianceReal :: f (a, Vector emiSh a) -> HermitianMatrix emiSh a} + +covarianceReal :: + (Shape.C emiSh, Eq emiSh, Class.Real a) => + NonEmpty.T [] (a, Vector emiSh a) -> HermitianMatrix emiSh a +covarianceReal = + getCovarianceReal $ + Class.switchReal + (CovarianceReal $ Hermitian.sumRank1NonEmpty MatrixShape.RowMajor) + (CovarianceReal $ Hermitian.sumRank1NonEmpty MatrixShape.RowMajor) + +gaussian :: + (Shape.C emiSh, Shape.C stateSh, Class.Real prob) => + Array stateSh (Vector emiSh prob, HermitianMatrix emiSh prob) -> + Gaussian emiSh stateSh prob +gaussian = consGaussian . fmap gaussianParameters + +gaussianParameters :: + (Shape.C emiSh, Class.Real prob) => + (Vector emiSh prob, HermitianMatrix emiSh prob) -> + (Vector emiSh prob, UpperTriangular emiSh prob, prob) +gaussianParameters (center, covariance) = + gaussianFromCholesky center $ HermitianPD.decompose covariance + +consGaussian :: + (Shape.C stateSh) => + Array stateSh (Vector emiSh a, UpperTriangular emiSh a, a) -> + Gaussian emiSh stateSh a +consGaussian = Gaussian + +gaussianFromCholesky :: + (Shape.C emiSh, Class.Real prob) => + Vector emiSh prob -> UpperTriangular emiSh prob -> + (Vector emiSh prob, UpperTriangular emiSh prob, prob) +gaussianFromCholesky center covarianceChol = + let covarianceSqrtDet = + Vector.product $ Triangular.takeDiagonal covarianceChol + in (center, covarianceChol, + recip (sqrt2pi ^ vectorDim center * covarianceSqrtDet)) + +sqrt2pi :: (Class.Real a) => a +sqrt2pi = runIdentity $ Class.switchReal sqrt2piAux sqrt2piAux + +sqrt2piAux :: (Floating a) => Identity a +sqrt2piAux = Identity $ sqrt (2*pi) + +cexp :: (Class.Real a) => a -> a +cexp = appEndo $ Class.switchReal (Endo exp) (Endo exp) + + + +class ToCSV distr where + toCells :: distr -> [[String]] + +class FromCSV distr where + parseCells :: StateShape distr -> HMMCSV.CSVParser distr + +class (Ord symbol) => CSVSymbol symbol where + cellFromSymbol :: symbol -> String + symbolFromCell :: String -> Maybe symbol + +instance CSVSymbol Char where + cellFromSymbol = (:[]) + symbolFromCell = listToMaybe + +instance CSVSymbol Int where + cellFromSymbol = show + symbolFromCell = maybeRead + + +instance + (Shape.C sh, Class.Real prob, Show prob, Read prob, CSVSymbol symbol) => + ToCSV (Discrete symbol sh prob) where + toCells (Discrete m) = + map + (\(symbol, probs) -> + cellFromSymbol symbol : HMMCSV.cellsFromVector probs) $ + Map.toAscList m + +instance + (Shape.C sh, Class.Real prob, Show prob, Read prob, CSVSymbol symbol) => + FromCSV (Discrete symbol sh prob) where + parseCells n = + fmap (Discrete . Map.fromList) $ + HMMCSV.manyRowsUntilEnd $ parseSymbolProb n + +parseSymbolProb :: + (Shape.C sh, Class.Real prob, Read prob, CSVSymbol symbol) => + sh -> CSV.CSVRow -> HMMCSV.CSVParser (symbol, Vector sh prob) +parseSymbolProb sh row = + case row of + [] -> MT.lift $ ME.throw "missing symbol" + c:cs -> + liftM2 (,) + (let str = CSV.csvFieldContent c + in MT.lift $ ME.fromMaybe (printf "unknown symbol %s" str) $ + symbolFromCell str) + (do v <- HMMCSV.parseVectorFields cs + let n = Shape.size sh + let m = vectorDim v + HMMCSV.assert (n == m) + (printf "number of states (%d) and size of probability vector (%d) mismatch" + n m) + return $ StorableArray.reshape sh v) + + +instance + (Shape.Indexed emiSh, Shape.Indexed stateSh, + Class.Real a, Eq a, Show a, Read a) => + ToCSV (Gaussian emiSh stateSh a) where + toCells (Gaussian params) = + List.intercalate [[]] $ + map + (\(center, covarianceChol, _) -> + HMMCSV.cellsFromVector center : + HMMCSV.cellsFromSquare (Triangular.toSquare covarianceChol)) $ + Array.toList params + +instance + (emiSh ~ Matrix.ZeroInt, Shape.Indexed stateSh, + Class.Real a, Eq a, Show a, Read a) => + FromCSV (Gaussian emiSh stateSh a) where + parseCells sh = do + let n = Shape.size sh + gs <- HMMCSV.manySepUntilEnd parseSingleGaussian + HMMCSV.assert (length gs == n) $ + printf "number of states (%d) and number of Gaussians (%d) mismatch" + n (length gs) + let sizes = map (vectorDim . fst3) gs + HMMCSV.assert (ListHT.allEqual sizes) $ + printf "dimensions of emissions mismatch: %s" (show sizes) + return $ consGaussian $ Array.fromList sh gs + +parseSingleGaussian :: + (emiSh ~ Matrix.ZeroInt, Class.Real prob, Eq prob, Read prob) => + HMMCSV.CSVParser (Vector emiSh prob, UpperTriangular emiSh prob, prob) +parseSingleGaussian = do + center <- HMMCSV.parseNonEmptyVectorCells + covarianceCholSquare <- + HMMCSV.parseSquareMatrixCells $ StorableArray.shape center + let covarianceChol = Triangular.takeUpper covarianceCholSquare + HMMCSV.assert + (isUpperTriang covarianceCholSquare covarianceChol) + "matrices must be upper triangular" + return $ gaussianFromCholesky center covarianceChol + + +{- +Maybe this test is too strict. +It would also be ok, and certainly more intuitive +to use an orthogonal but not normalized matrix. +We could get such a matrix from the eigensystem. +-} +isUpperTriang :: + (Shape.C sh, Class.Real a, Eq a) => + SquareMatrix sh a -> UpperTriangular sh a -> Bool +isUpperTriang m mt = + Vector.toList m == Vector.toList (Triangular.toSquare mt) diff --git a/src/Math/HiddenMarkovModel/Example/Circle.hs b/src/Math/HiddenMarkovModel/Example/Circle.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Example/Circle.hs @@ -0,0 +1,12 @@ +{- | +Example of an HMM with continuous emissions with two-dimensional observations. +We train a model to accept a parametric curve of a circle with a certain speed. +This is like "Math.HiddenMarkovModel.Example.SineWave" but in two dimensions. + +The four hidden states correspond to the four quadrants. +-} +module Math.HiddenMarkovModel.Example.Circle +{-# WARNING "do not import that module, it is only intended for demonstration" #-} + (module Math.HiddenMarkovModel.Example.CirclePrivate) where + +import Math.HiddenMarkovModel.Example.CirclePrivate diff --git a/src/Math/HiddenMarkovModel/Example/CirclePrivate.hs b/src/Math/HiddenMarkovModel/Example/CirclePrivate.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Example/CirclePrivate.hs @@ -0,0 +1,123 @@ +module Math.HiddenMarkovModel.Example.CirclePrivate where + +import qualified Math.HiddenMarkovModel as HMM +import qualified Math.HiddenMarkovModel.Distribution as Distr +import Math.HiddenMarkovModel.Utility + (normalizeProb, squareFromLists, hermitianFromList) + +import qualified Numeric.LAPACK.Vector as Vector +import Numeric.LAPACK.Vector (Vector) + +import qualified Data.Array.Comfort.Boxed as Array +import qualified Data.Array.Comfort.Shape as Shape + +import qualified System.Random as Rnd + +import qualified Control.Monad.Trans.State as MS +import Control.Monad (liftM2, replicateM) + +import qualified Data.NonEmpty.Class as NonEmptyC +import qualified Data.NonEmpty as NonEmpty +import Data.Function.HT (nest) +import Data.NonEmpty ((!:)) +import Data.Maybe (fromMaybe) + + + +data State = Q1 | Q2 | Q3 | Q4 + deriving (Eq, Ord, Enum, Bounded) + +type StateSet = Shape.Enumeration State + +stateSet :: StateSet +stateSet = Shape.Enumeration + + +data Coordinate = X | Y + deriving (Eq, Ord, Enum, Bounded) + +type CoordinateSet = Shape.Enumeration Coordinate + +coordinateSet :: CoordinateSet +coordinateSet = Shape.Enumeration + +type HMM = HMM.Gaussian CoordinateSet StateSet Double + +hmm :: HMM +hmm = + HMM.Cons { + HMM.initial = normalizeProb $ Vector.constant stateSet 1, + HMM.transition = + squareFromLists stateSet $ + stateVector 0.9 0.0 0.0 0.1 : + stateVector 0.1 0.9 0.0 0.0 : + stateVector 0.0 0.1 0.9 0.0 : + stateVector 0.0 0.0 0.1 0.9 : + [], + HMM.distribution = + let cov0 = hermitianFromList coordinateSet [0.10, -0.09, 0.10] + cov1 = hermitianFromList coordinateSet [0.10, 0.09, 0.10] + in Distr.gaussian $ Array.fromList stateSet $ + (Vector.fromList coordinateSet [ 0.5, 0.5], cov0) : + (Vector.fromList coordinateSet [-0.5, 0.5], cov1) : + (Vector.fromList coordinateSet [-0.5, -0.5], cov0) : + (Vector.fromList coordinateSet [ 0.5, -0.5], cov1) : + [] + } + +stateVector :: Double -> Double -> Double -> Double -> Vector StateSet Double +stateVector x0 x1 x2 x3 = Vector.fromList stateSet [x0,x1,x2,x3] + +circleLabeled :: NonEmpty.T [] (State, Vector CoordinateSet Double) +circleLabeled = + NonEmpty.mapTail (take 200) $ + fmap + (\x -> + (toEnum $ mod (floor (x*2/pi)) 4, + Vector.fromList coordinateSet [cos x, sin x])) $ + NonEmptyC.iterate (0.1+) 0 + +circle :: NonEmpty.T [] (Vector CoordinateSet Double) +circle = fmap snd circleLabeled + +revealed :: NonEmpty.T [] State +revealed = HMM.reveal hmm circle + +{- | +Sample multivariate normal distribution and reconstruct it from the samples. +You should obtain the same parameters. +-} +reconstructDistribution :: HMM.Gaussian CoordinateSet () Double +reconstructDistribution = + let gen = Distr.generate (HMM.distribution hmm) Q1 + in HMM.finishTraining $ HMM.trainSupervised () $ fmap ((,) ()) $ + flip MS.evalState (Rnd.mkStdGen 23) $ + liftM2 (!:) gen $ replicateM 1000 gen + +{- | +Generate labeled emission sequences +and use them for supervised training. +-} +reconstructModel :: HMM +reconstructModel = + HMM.trainMany (HMM.trainSupervised stateSet) $ + fmap + (\seed -> + fromMaybe (error "empty generated sequence") $ NonEmpty.fetch $ + take 1000 $ HMM.generateLabeled hmm $ Rnd.mkStdGen seed) + (23 !: take 42 [24..]) + + +hmmTrainedSupervised :: HMM +hmmTrainedSupervised = + HMM.finishTraining $ HMM.trainSupervised stateSet circleLabeled + +hmmTrainedUnsupervised :: HMM +hmmTrainedUnsupervised = + HMM.finishTraining $ HMM.trainUnsupervised hmm circle + +hmmIterativelyTrained :: HMM +hmmIterativelyTrained = + nest 100 + (HMM.finishTraining . flip HMM.trainUnsupervised circle) + hmm diff --git a/src/Math/HiddenMarkovModel/Example/SineWave.hs b/src/Math/HiddenMarkovModel/Example/SineWave.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Example/SineWave.hs @@ -0,0 +1,90 @@ +{- | +Example of an HMM with continuous emissions. +We train a model to accept sine waves of a certain frequency. + +There are four hidden states: 'Rising', 'High', 'Falling', 'Low'. +-} +module Math.HiddenMarkovModel.Example.SineWave +{-# WARNING "do not import that module, it is only intended for demonstration" #-} + where + +import qualified Math.HiddenMarkovModel as HMM +import qualified Math.HiddenMarkovModel.Distribution as Distr +import Math.HiddenMarkovModel.Utility + (normalizeProb, squareFromLists, hermitianFromList, singleton) + +import qualified Numeric.LAPACK.Vector as Vector +import Numeric.LAPACK.Vector (Vector) + +import qualified Data.Array.Comfort.Boxed as Array +import qualified Data.Array.Comfort.Shape as Shape + +import qualified Data.NonEmpty.Class as NonEmptyC +import qualified Data.NonEmpty as NonEmpty +import Data.Function.HT (nest) +import Data.Tuple.HT (mapSnd) + + + +data State = Rising | High | Falling | Low + deriving (Eq, Ord, Enum, Bounded) + +type StateSet = Shape.Enumeration State + +stateSet :: StateSet +stateSet = Shape.Enumeration + + +type HMM = HMM.Gaussian () StateSet Double + +hmm :: HMM +hmm = + HMM.Cons { + HMM.initial = normalizeProb $ Vector.constant stateSet 1, + HMM.transition = + squareFromLists stateSet $ + stateVector 0.9 0.0 0.0 0.1 : + stateVector 0.1 0.9 0.0 0.0 : + stateVector 0.0 0.1 0.9 0.0 : + stateVector 0.0 0.0 0.1 0.9 : + [], + HMM.distribution = + Distr.gaussian $ Array.fromList stateSet $ + (singleton 0 , hermitianFromList () [1]) : + (singleton 1 , hermitianFromList () [1]) : + (singleton 0 , hermitianFromList () [1]) : + (singleton (-1), hermitianFromList () [1]) : + [] + } + +stateVector :: Double -> Double -> Double -> Double -> Vector StateSet Double +stateVector x0 x1 x2 x3 = Vector.fromList stateSet [x0,x1,x2,x3] + +sineWaveLabeled :: NonEmpty.T [] (State, Double) +sineWaveLabeled = + NonEmpty.mapTail (take 200) $ + fmap (\x -> (toEnum $ mod (floor (x*2/pi+0.5)) 4, sin x)) $ + NonEmptyC.iterate (0.1+) 0 + +sineWave :: NonEmpty.T [] Double +sineWave = fmap snd sineWaveLabeled + +revealed :: NonEmpty.T [] State +revealed = HMM.reveal hmmTrainedSupervised $ fmap singleton sineWave + +hmmTrainedSupervised :: HMM +hmmTrainedSupervised = + HMM.finishTraining $ HMM.trainSupervised stateSet $ + fmap (mapSnd singleton) sineWaveLabeled + +hmmTrainedUnsupervised :: HMM +hmmTrainedUnsupervised = + HMM.finishTraining $ HMM.trainUnsupervised hmm $ fmap singleton sineWave + +hmmIterativelyTrained :: HMM +hmmIterativelyTrained = + nest 100 + (\model -> + HMM.finishTraining $ HMM.trainUnsupervised model $ + fmap singleton sineWave) + hmm diff --git a/src/Math/HiddenMarkovModel/Example/TrafficLight.hs b/src/Math/HiddenMarkovModel/Example/TrafficLight.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Example/TrafficLight.hs @@ -0,0 +1,50 @@ +{- | +This is an example of an HMM with discrete emissions. +We model a traffic light consisting of the colors red, yellow, green, +where only one lamp can be switched on at every point in time. +This way, when it is yellow you cannot tell immediately +whether it will switch to green or red. +We can only infer this from the light seen before. + +There are four hidden states: +'StateRed' emits red, 'StateYellowRG' emits yellow between red and green, +'StateGreen' emits green, 'StateYellowGR' emits yellow between green and red. + +We quantise time in time steps. +The transition matrix of the model 'hmm' encodes +the expected duration of every state counted in time steps +and what states follow after each other. +E.g. transition probability of 0.8 of a state to itself means +that the expected duration of the state is 5 time steps (1/(1-0.8)). +However, it is a geometric distribution, +that is, shorter durations are always more probable. + +The distribution of 'hmm' encodes which lights a state activates. +In our case everything is deterministic: +Every state can switch exactly one light on. + +Given a sequence of observed lights +the function 'HMM.reveal' tells us the most likely sequence of states. +We test this with the light sequences in 'stateSequences' +where we already know the hidden states +as they are stored in 'labeledSequences'. +'verifyRevelation' compares the computed state sequence with the given one. + +We also try some trainings in 'hmmTrainedSupervised' et.al. +-} +module Math.HiddenMarkovModel.Example.TrafficLight +{-# WARNING "do not import that module, it is only intended for demonstration" #-} + ( + HMM, + Color(..), + hmm, + hmmDisturbed, + red, yellowRG, green, yellowGR, + labeledSequences, + hmmTrainedSupervised, + stateSequences, + hmmTrainedUnsupervised, + hmmIterativelyTrained, + ) where + +import Math.HiddenMarkovModel.Example.TrafficLightPrivate diff --git a/src/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs b/src/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs @@ -0,0 +1,164 @@ +{-# LANGUAGE TypeFamilies #-} +module Math.HiddenMarkovModel.Example.TrafficLightPrivate where + +import qualified Math.HiddenMarkovModel as HMM +import qualified Math.HiddenMarkovModel.Distribution as Distr +import Math.HiddenMarkovModel.Utility (normalizeProb, squareFromLists) + +import qualified Numeric.LAPACK.Vector as Vector +import Numeric.LAPACK.Vector (Vector) +import Numeric.LAPACK.Format (Format(format)) + +import qualified Data.Array.Comfort.Shape as Shape + +import qualified Text.PrettyPrint.Boxes as TextBox +import Text.Read.HT (maybeRead) + +import Control.DeepSeq (NFData(rnf)) +import Control.Monad (liftM2) + +import qualified Data.Map as Map +import qualified Data.NonEmpty as NonEmpty +import qualified Data.List.HT as ListHT +import Data.NonEmpty ((!:)) + + + +data Color = Red | Yellow | Green + deriving (Eq, Ord, Enum, Show, Read) + +instance NFData Color where + rnf Red = () + rnf _ = () + +instance Format Color where + format _fmt = TextBox.text . show + +{- | +Using 'show' and 'read' is not always a good choice +since they must format and parse Haskell expressions +which is not of much use to the outside world. +-} +instance Distr.CSVSymbol Color where + cellFromSymbol = show + symbolFromCell = maybeRead + + +-- data State = StateRed | StateYellowDown | StateGreen | StateYellowUp +data State = StateRed | StateYellowRG | StateGreen | StateYellowGR + deriving (Eq, Ord, Enum, Bounded) + +type StateSet = Shape.Enumeration State + +stateSet :: StateSet +stateSet = Shape.Enumeration + + +type HMM = HMM.Discrete Color StateSet Double + +hmm :: HMM +hmm = + HMM.Cons { + HMM.initial = normalizeProb $ stateVector 2 1 2 1, + HMM.transition = + squareFromLists stateSet $ + stateVector 0.8 0.0 0.0 0.2 : + stateVector 0.2 0.8 0.0 0.0 : + stateVector 0.0 0.2 0.8 0.0 : + stateVector 0.0 0.0 0.2 0.8 : + [], + HMM.distribution = + Distr.Discrete $ Map.fromList $ + (Red, stateVector 1 0 0 0) : + (Yellow, stateVector 0 1 0 1): + (Green, stateVector 0 0 1 0) : + [] + } + +hmmDisturbed :: HMM +hmmDisturbed = + HMM.Cons { + HMM.initial = normalizeProb $ stateVector 1 1 1 1, + HMM.transition = + squareFromLists stateSet $ + stateVector 0.3 0.2 0.2 0.3 : + stateVector 0.3 0.3 0.2 0.2 : + stateVector 0.2 0.3 0.3 0.2 : + stateVector 0.2 0.2 0.3 0.3 : + [], + HMM.distribution = + Distr.Discrete $ Map.fromList $ + (Red, stateVector 0.6 0.2 0.2 0.2) : + (Yellow, stateVector 0.2 0.6 0.2 0.6) : + (Green, stateVector 0.2 0.2 0.6 0.2) : + [] + } + +stateVector :: Double -> Double -> Double -> Double -> Vector StateSet Double +stateVector x0 x1 x2 x3 = Vector.fromList stateSet [x0,x1,x2,x3] + + +red, yellowRG, green, yellowGR :: (State, Color) +red = (StateRed, Red) +yellowRG = (StateYellowRG, Yellow) +green = (StateGreen, Green) +yellowGR = (StateYellowGR, Yellow) + +labeledSequences :: NonEmpty.T [] (NonEmpty.T [] (State, Color)) +labeledSequences = + (red !: red : red : red : + yellowRG : yellowRG : + green : green : green : green : green : + yellowGR : + red : red : red : + []) !: + (green !: green : green : + yellowGR : + red : red : red : red : + yellowRG : + green : green : green : green : green : + yellowGR : yellowGR : + []) : + [] + +{- | +Construct a Hidden Markov model by watching a set +of manually created sequences of emissions and according states. +-} +hmmTrainedSupervised :: HMM +hmmTrainedSupervised = + HMM.trainMany (HMM.trainSupervised stateSet) labeledSequences + + +stateSequences :: NonEmpty.T [] (NonEmpty.T [] Color) +stateSequences = fmap (fmap snd) labeledSequences + +{- | +Construct a Hidden Markov model starting from a known model +and a set of sequences that contain only the emissions, but no states. +-} +hmmTrainedUnsupervised :: HMM +hmmTrainedUnsupervised = + HMM.trainMany (HMM.trainUnsupervised hmm) stateSequences + +{- | +Repeat unsupervised training until convergence. +-} +hmmIterativelyTrained :: HMM +hmmIterativelyTrained = + snd $ head $ dropWhile fst $ + ListHT.mapAdjacent (\hmm0 hmm1 -> (HMM.deviation hmm0 hmm1 > 1e-5, hmm1)) $ + iterate + (flip HMM.trainMany stateSequences . HMM.trainUnsupervised) + hmmDisturbed + + +verifyRevelation :: HMM -> NonEmpty.T [] (State, Color) -> Bool +verifyRevelation model xs = + fmap fst xs == HMM.reveal model (fmap snd xs) + +verifyRevelations :: [Bool] +verifyRevelations = + liftM2 verifyRevelation + [hmm, hmmDisturbed, hmmTrainedSupervised, hmmTrainedUnsupervised] + (NonEmpty.flatten labeledSequences) diff --git a/src/Math/HiddenMarkovModel/Named.hs b/src/Math/HiddenMarkovModel/Named.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Named.hs @@ -0,0 +1,117 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE UndecidableInstances #-} +module Math.HiddenMarkovModel.Named ( + T(..), + Discrete, + Gaussian, + fromModelAndNames, + toCSV, + fromCSV, + ) where + +import qualified Math.HiddenMarkovModel.Distribution as Distr +import qualified Math.HiddenMarkovModel.Private as HMM +import qualified Math.HiddenMarkovModel.CSV as HMMCSV +import Math.HiddenMarkovModel.Utility (attachOnes, vectorDim) + +import qualified Numeric.Netlib.Class as Class + +import qualified Data.Array.Comfort.Storable as StorableArray +import qualified Data.Array.Comfort.Boxed as Array +import qualified Data.Array.Comfort.Shape as Shape +import Data.Array.Comfort.Boxed (Array) + +import qualified Text.CSV.Lazy.String as CSV +import Text.Printf (printf) + +import qualified Control.Monad.Exception.Synchronous as ME +import qualified Control.Monad.Trans.State as MS +import Control.DeepSeq (NFData, rnf) +import Foreign.Storable (Storable) + +import qualified Data.Map as Map +import qualified Data.List as List +import Data.Tuple.HT (swap) +import Data.Map (Map) + + +{- | +A Hidden Markov Model with names for each state. + +Although 'nameFromStateMap' and 'stateFromNameMap' are exported +you must be careful to keep them consistent when you alter them. +-} +data T distr sh ix prob = + Cons { + model :: HMM.T distr sh prob, + nameFromStateMap :: Array sh String, + stateFromNameMap :: Map String ix + } + deriving (Show) + +type Discrete symbol stateSh prob = + T (Distr.Discrete symbol stateSh prob) stateSh (Shape.Index stateSh) prob +type Gaussian emiSh stateSh a = + T (Distr.Gaussian emiSh stateSh a) stateSh (Shape.Index stateSh) a + + +instance + (NFData distr, NFData sh, NFData ix, NFData prob, + Shape.C sh, Storable prob) => + NFData (T distr sh ix prob) where + rnf hmm = rnf (model hmm, nameFromStateMap hmm, stateFromNameMap hmm) + + +fromModelAndNames :: + (Shape.Indexed sh, Shape.Index sh ~ state) => + HMM.T distr sh prob -> [String] -> T distr sh state prob +fromModelAndNames md names = + let m = Array.fromList (StorableArray.shape $ HMM.initial md) names + in Cons { + model = md, + nameFromStateMap = m, + stateFromNameMap = inverseMap m + } + +inverseMap :: + (Shape.Indexed sh, Shape.Index sh ~ ix) => Array sh String -> Map String ix +inverseMap = + Map.fromListWith (error "duplicate label") . + map swap . Array.toAssociations + + +toCSV :: + (Distr.ToCSV distr, Shape.Indexed sh, Class.Real prob, Show prob) => + T distr sh ix prob -> String +toCSV hmm = + CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $ + Array.toList (nameFromStateMap hmm) : HMM.toCells (model hmm) + +fromCSV :: + (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh, + Shape.Indexed stateSh, Shape.Index stateSh ~ state, + Class.Real prob, Read prob) => + (Int -> stateSh) -> + String -> ME.Exceptional String (T distr stateSh state prob) +fromCSV makeShape = + MS.evalStateT (parseCSV makeShape) . map HMMCSV.fixShortRow . CSV.parseCSV + +parseCSV :: + (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh, + Shape.Indexed stateSh, Shape.Index stateSh ~ state, + Class.Real prob, Read prob) => + (Int -> stateSh) -> HMMCSV.CSVParser (T distr stateSh state prob) +parseCSV makeShape = do + names <- HMMCSV.parseStringList =<< HMMCSV.getRow + let duplicateNames = + Map.keys $ Map.filter (> (1::Int)) $ + Map.fromListWith (+) $ attachOnes names + in HMMCSV.assert (null duplicateNames) $ + "duplicate names: " ++ List.intercalate ", " duplicateNames + md <- HMM.parseCSV makeShape + let n = length names + m = vectorDim (HMM.initial md) + in HMMCSV.assert (n == m) $ + printf "got %d state names for %d states" n m + return $ fromModelAndNames md names diff --git a/src/Math/HiddenMarkovModel/Normalized.hs b/src/Math/HiddenMarkovModel/Normalized.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Normalized.hs @@ -0,0 +1,175 @@ +{-# LANGUAGE TypeFamilies #-} +{- | +Counterparts to functions in "Math.HiddenMarkovModel.Private" +that normalize interim results. +We need to do this in order to prevent +to round very small probabilities to zero. +-} +module Math.HiddenMarkovModel.Normalized where + +import qualified Math.HiddenMarkovModel.Distribution as Distr +import Math.HiddenMarkovModel.Private + (T(..), Trained(..), emission, + biscaleTransition, matrixMaxMul, sumTransitions) +import Math.HiddenMarkovModel.Utility + (SquareMatrix, normalizeFactor, normalizeProb) + +import qualified Numeric.LAPACK.Vector as Vector +import Numeric.LAPACK.Matrix ((<#), (#>)) +import Numeric.LAPACK.Vector (Vector) + +import qualified Numeric.Netlib.Class as Class + +import qualified Control.Functor.HT as Functor + +import qualified Data.Array.Comfort.Storable as StorableArray +import qualified Data.Array.Comfort.Boxed as Array +import qualified Data.Array.Comfort.Shape as Shape + +import qualified Data.NonEmpty.Class as NonEmptyC +import qualified Data.NonEmpty as NonEmpty +import qualified Data.Foldable as Fold +import qualified Data.List as List +import Data.Traversable (Traversable, mapAccumL) +import Data.Tuple.HT (mapFst, mapSnd, swap) + + +{- | +Logarithm of the likelihood to observe the given sequence. +We return the logarithm because the likelihood can be so small +that it may be rounded to zero in the choosen number type. +-} +logLikelihood :: + (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh, Floating prob, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr sh prob -> NonEmpty.T f emission -> prob +logLikelihood hmm = Fold.sum . fmap (log . fst) . alpha hmm + +alpha :: + (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr sh prob -> + NonEmpty.T f emission -> NonEmpty.T f (prob, Vector sh prob) +alpha hmm (NonEmpty.Cons x xs) = + let normMulEmiss y = normalizeFactor . Vector.mul (emission hmm y) + in NonEmpty.scanl + (\(_,alphai) xi -> normMulEmiss xi (transition hmm #> alphai)) + (normMulEmiss x (initial hmm)) + xs + +beta :: + (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f, NonEmptyC.Reverse f) => + T distr sh prob -> + f (prob, emission) -> NonEmpty.T f (Vector sh prob) +beta hmm = + nonEmptyScanr + (\(ci,xi) betai -> + Vector.scale (recip ci) $ + Vector.mul (emission hmm xi) betai <# transition hmm) + (Vector.constant (StorableArray.shape $ initial hmm) 1) + +alphaBeta :: + (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f, NonEmptyC.Zip f, NonEmptyC.Reverse f) => + T distr sh prob -> + NonEmpty.T f emission -> + (NonEmpty.T f (prob, Vector sh prob), NonEmpty.T f (Vector sh prob)) +alphaBeta hmm xs = + let calphas = alpha hmm xs + in (calphas, + beta hmm $ NonEmpty.tail $ NonEmptyC.zip (fmap fst calphas) xs) + + +xiFromAlphaBeta :: + (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f, NonEmptyC.Zip f) => + T distr sh prob -> + NonEmpty.T f emission -> + NonEmpty.T f (prob, Vector sh prob) -> + NonEmpty.T f (Vector sh prob) -> + f (SquareMatrix sh prob) +xiFromAlphaBeta hmm xs calphas betas = + let (cs,alphas) = Functor.unzip calphas + in NonEmptyC.zipWith4 + (\x alpha0 c1 beta1 -> + Vector.scale (recip c1) $ biscaleTransition hmm x alpha0 beta1) + (NonEmpty.tail xs) + (NonEmpty.init alphas) + (NonEmpty.tail cs) + (NonEmpty.tail betas) + +zetaFromAlphaBeta :: + (Shape.C sh, Eq sh, Class.Real prob, NonEmptyC.Zip f) => + NonEmpty.T f (prob, Vector sh prob) -> + NonEmpty.T f (Vector sh prob) -> + NonEmpty.T f (Vector sh prob) +zetaFromAlphaBeta calphas betas = + NonEmptyC.zipWith (Vector.mul . snd) calphas betas + + +{- | +Reveal the state sequence +that led most likely to the observed sequence of emissions. +It is found using the Viterbi algorithm. +-} +reveal :: + (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, + Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ state, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f, NonEmptyC.Reverse f) => + T distr sh prob -> NonEmpty.T f emission -> NonEmpty.T f state +reveal hmm (NonEmpty.Cons x xs) = + fmap (Shape.revealIndex (StorableArray.shape $ initial hmm)) $ + uncurry (NonEmpty.scanr (StorableArray.!)) $ + mapFst + (fst . Vector.argAbsMaximum . + StorableArray.mapShape Shape.Deferred) $ + mapAccumL + (\alphai xi -> + swap $ mapSnd (Vector.mul (emission hmm xi)) $ + matrixMaxMul (transition hmm) $ normalizeProb alphai) + (Vector.mul (emission hmm x) (initial hmm)) xs + + +{- | +Variant of NonEmpty.scanr with less stack consumption. +-} +nonEmptyScanr :: + (Traversable f, NonEmptyC.Reverse f) => + (a -> b -> b) -> b -> f a -> NonEmpty.T f b +nonEmptyScanr f x = + NonEmptyC.reverse . NonEmpty.scanl (flip f) x . NonEmptyC.reverse + + +{- | +Consider a superposition of all possible state sequences +weighted by the likelihood to produce the observed emission sequence. +Now train the model with respect to all of these sequences +with respect to the weights. +This is done by the Baum-Welch algorithm. +-} +trainUnsupervised :: + (Distr.Estimate tdistr distr, Distr.StateShape distr ~ sh, Eq sh, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + T distr sh prob -> NonEmpty.T [] emission -> Trained tdistr sh prob +trainUnsupervised hmm xs = + let (alphas, betas) = alphaBeta hmm xs + zetas = zetaFromAlphaBeta alphas betas + zeta0 = NonEmpty.head zetas + + in Trained { + trainedInitial = zeta0, + trainedTransition = + sumTransitions hmm $ xiFromAlphaBeta hmm xs alphas betas, + trainedDistribution = + Distr.accumulateEmissions $ + Array.fromList (StorableArray.shape zeta0) $ + map (zip (NonEmpty.flatten xs)) $ + List.transpose $ map Vector.toList $ NonEmpty.flatten zetas + } diff --git a/src/Math/HiddenMarkovModel/Pattern.hs b/src/Math/HiddenMarkovModel/Pattern.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Pattern.hs @@ -0,0 +1,121 @@ +{-# LANGUAGE TypeFamilies #-} +{- | +This module provides a simple way to train +the transition matrix and initial probability vector +using simple patterns of state sequences. + +You may create a trained model using semigroup combinators like this: + +> example :: HMM.DiscreteTrained Char (ShapeStatic.ZeroBased TypeNum.U2) Double +> example = +> let a = atom FL.i0 +> b = atom FL.i1 +> distr = +> Distr.DiscreteTrained $ Map.fromList $ +> ('a', ShapeStatic.vector $ 1!:2!:FL.end) : +> ('b', ShapeStatic.vector $ 4!:3!:FL.end) : +> ('c', ShapeStatic.vector $ 0!:1!:FL.end) : +> [] +> in finish (ShapeStatic.ZeroBased Proxy) distr $ +> replicate 5 $ replicate 10 a <> replicate 20 b +-} +module Math.HiddenMarkovModel.Pattern ( + T, + atom, + append, + replicate, + finish, + ) where + +import qualified Math.HiddenMarkovModel.Distribution as Distr +import qualified Math.HiddenMarkovModel as HMM +import Math.HiddenMarkovModel.Private (Trained(..)) +import Math.HiddenMarkovModel.Utility (SquareMatrix, squareConstant) + +import qualified Numeric.LAPACK.Vector as Vector +import qualified Numeric.LAPACK.ShapeStatic as ShapeStatic + +import qualified Numeric.Netlib.Class as Class + +import qualified Data.Array.Comfort.Storable as StorableArray +import qualified Data.Array.Comfort.Shape as Shape + +import qualified Data.FixedLength as FL +import Data.FixedLength ((!:)) + +import qualified Type.Data.Num.Unary.Literal as TypeNum +import Type.Base.Proxy (Proxy(Proxy)) + +import qualified Data.Map as Map +import Data.Semigroup (Semigroup, (<>), stimes) + +import Prelude hiding (replicate) + + +newtype T sh prob = + Cons (sh -> (Shape.Index sh, SquareMatrix sh prob, Shape.Index sh)) + +atom :: + (Shape.Indexed sh, Shape.Index sh ~ state, Class.Real prob) => + state -> T sh prob +atom s = Cons $ \sh -> (s, squareConstant sh 0, s) + + +instance + (Shape.Indexed sh, Eq sh, Class.Real prob) => + Semigroup (T sh prob) where + (<>) = append + stimes k = replicate $ fromIntegral k + + +infixl 5 `append` + +append :: + (Shape.Indexed sh, Eq sh, Class.Real prob) => + T sh prob -> T sh prob -> T sh prob +append (Cons f) (Cons g) = + Cons $ \n -> + case (f n, g n) of + ((sai, ma, sao), (sbi, mb, sbo)) -> + (sai, increment (sbi,sao) 1 $ Vector.add ma mb, sbo) + +replicate :: + (Shape.Indexed sh, Class.Real prob) => Int -> T sh prob -> T sh prob +replicate ki (Cons f) = + Cons $ \sh -> + case f sh of + (si, m, so) -> + let k = fromIntegral ki + in (si, increment (si,so) (k-1) $ Vector.scale k m, so) + +increment :: + (Shape.Indexed sh, Shape.Index sh ~ state, Class.Real a) => + (state, state) -> a -> SquareMatrix sh a -> SquareMatrix sh a +increment (i,j) x m = StorableArray.accumulate (+) m [((i,j), x)] + + +finish :: + (Shape.Indexed sh, Class.Real prob) => + sh -> tdistr -> T sh prob -> Trained tdistr sh prob +finish sh tdistr (Cons f) = + case f sh of + (si, m, _so) -> + Trained { + trainedInitial = StorableArray.fromAssociations sh 0 [(si,1)], + trainedTransition = m, + trainedDistribution = tdistr + } + + +_example :: HMM.DiscreteTrained Char (ShapeStatic.ZeroBased TypeNum.U2) Double +_example = + let a = atom FL.i0 + b = atom FL.i1 + distr = + Distr.DiscreteTrained $ Map.fromList $ + ('a', ShapeStatic.vector $ 1!:2!:FL.end) : + ('b', ShapeStatic.vector $ 4!:3!:FL.end) : + ('c', ShapeStatic.vector $ 0!:1!:FL.end) : + [] + in finish (ShapeStatic.ZeroBased Proxy) distr $ + replicate 5 $ replicate 10 a <> replicate 20 b diff --git a/src/Math/HiddenMarkovModel/Private.hs b/src/Math/HiddenMarkovModel/Private.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Private.hs @@ -0,0 +1,335 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UndecidableInstances #-} +module Math.HiddenMarkovModel.Private where + +import qualified Math.HiddenMarkovModel.Distribution as Distr +import qualified Math.HiddenMarkovModel.CSV as HMMCSV +import Math.HiddenMarkovModel.Utility (SquareMatrix, diagonal) + +import qualified Numeric.LAPACK.Matrix as Matrix +import qualified Numeric.LAPACK.Vector as Vector +import qualified Numeric.LAPACK.Format as Format +import Numeric.LAPACK.Matrix ((<#), (<#>), (#>)) +import Numeric.LAPACK.Vector (Vector) + +import qualified Numeric.Netlib.Class as Class + +import Control.DeepSeq (NFData, rnf) +import Control.Applicative ((<$>)) + +import Foreign.Storable (Storable) + +import qualified Data.Array.Comfort.Storable as StorableArray +import qualified Data.Array.Comfort.Boxed as Array +import qualified Data.Array.Comfort.Shape as Shape + +import qualified Data.NonEmpty.Class as NonEmptyC +import qualified Data.NonEmpty as NonEmpty +import qualified Data.Semigroup as Sg +import qualified Data.List as List +import Data.Traversable (Traversable, mapAccumL) +import Data.Tuple.HT (mapPair, mapFst, mapSnd, swap) + + +{- | +A Hidden Markov model consists of a number of (hidden) states +and a set of emissions. +There is a vector for the initial probability of each state +and a matrix containing the probability for switching +from one state to another one. +The 'distribution' field points to probability distributions +that associate every state with emissions of different probability. +Famous distribution instances are discrete and Gaussian distributions. +See "Math.HiddenMarkovModel.Distribution" for details. + +The transition matrix is transposed +with respect to popular HMM descriptions. +But I think this is the natural orientation, because this way +you can write \"transition matrix times probability column vector\". + +The type has two type parameters, +although the one for the distribution would be enough. +However, replacing @prob@ by @Distr.Probability distr@ +would prohibit the derived Show and Read instances. +-} +data T distr sh prob = + Cons { + initial :: Vector sh prob, + transition :: SquareMatrix sh prob, + distribution :: distr + } + deriving (Show) + +instance + (NFData distr, NFData sh, NFData prob, Storable prob) => + NFData (T distr sh prob) where + rnf (Cons initial_ transition_ distribution_) = + rnf (initial_, transition_, distribution_) + +instance + (Class.Real prob, Format.FormatArray sh, Format.Format distr) => + Format.Format (T distr sh prob) where + format fmt (Cons initial_ transition_ distribution_) = + Format.format fmt (initial_, transition_, distribution_) + + +emission :: + (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + T distr sh prob -> emission -> Vector sh prob +emission = Distr.emissionProb . distribution + + +forward :: + (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr sh prob -> NonEmpty.T f emission -> prob +forward hmm = Vector.sum . NonEmpty.last . alpha hmm + +alpha :: + (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr sh prob -> + NonEmpty.T f emission -> NonEmpty.T f (Vector sh prob) +alpha hmm (NonEmpty.Cons x xs) = + NonEmpty.scanl + (\alphai xi -> Vector.mul (emission hmm xi) (transition hmm #> alphai)) + (Vector.mul (emission hmm x) (initial hmm)) + xs + + +backward :: + (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr sh prob -> NonEmpty.T f emission -> prob +backward hmm (NonEmpty.Cons x xs) = + Vector.sum $ + Vector.mul (initial hmm) $ + Vector.mul (emission hmm x) $ + NonEmpty.head $ beta hmm xs + +beta :: + (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr sh prob -> + f emission -> NonEmpty.T f (Vector sh prob) +beta hmm = + NonEmpty.scanr + (\xi betai -> Vector.mul (emission hmm xi) betai <# transition hmm) + (Vector.constant (StorableArray.shape $ initial hmm) 1) + + +alphaBeta :: + (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr sh prob -> + NonEmpty.T f emission -> + (prob, NonEmpty.T f (Vector sh prob), NonEmpty.T f (Vector sh prob)) +alphaBeta hmm xs = + let alphas = alpha hmm xs + betas = beta hmm $ NonEmpty.tail xs + recipLikelihood = recip $ Vector.sum $ NonEmpty.last alphas + in (recipLikelihood, alphas, betas) + + + +biscaleTransition :: + (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, + Distr.EmissionProb distr, Distr.Probability distr ~ prob) => + T distr sh prob -> Distr.Emission distr -> + Vector sh prob -> Vector sh prob -> SquareMatrix sh prob +biscaleTransition hmm x alpha0 beta1 = + diagonal (Vector.mul (emission hmm x) beta1) + <#> + transition hmm + <#> + diagonal alpha0 + +xiFromAlphaBeta :: + (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + T distr sh prob -> prob -> + NonEmpty.T [] emission -> + NonEmpty.T [] (Vector sh prob) -> + NonEmpty.T [] (Vector sh prob) -> + [SquareMatrix sh prob] +xiFromAlphaBeta hmm recipLikelihood xs alphas betas = + zipWith3 + (\x alpha0 beta1 -> + Vector.scale recipLikelihood $ + biscaleTransition hmm x alpha0 beta1) + (NonEmpty.tail xs) + (NonEmpty.init alphas) + (NonEmpty.tail betas) + +zetaFromXi :: + (Shape.C sh, Eq sh, Class.Real prob) => + [SquareMatrix sh prob] -> [Vector sh prob] +zetaFromXi = map Matrix.columnSums + +zetaFromAlphaBeta :: + (Shape.C sh, Eq sh, Class.Real prob) => + prob -> + NonEmpty.T [] (Vector sh prob) -> + NonEmpty.T [] (Vector sh prob) -> + NonEmpty.T [] (Vector sh prob) +zetaFromAlphaBeta recipLikelihood alphas betas = + fmap (Vector.scale recipLikelihood) $ + NonEmptyC.zipWith Vector.mul alphas betas + + +{- | +In constrast to Math.HiddenMarkovModel.reveal +this does not normalize the vector. +This is slightly simpler but for long sequences +the product of probabilities might be smaller +than the smallest representable number. +-} +reveal :: + (Shape.InvIndexed sh, Shape.Index sh ~ state, + Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr sh prob -> NonEmpty.T f emission -> NonEmpty.T f state +reveal hmm (NonEmpty.Cons x xs) = + fmap (Shape.revealIndex (StorableArray.shape $ initial hmm)) $ + uncurry (NonEmpty.scanr (StorableArray.!)) $ + mapFst + (fst . Vector.argAbsMaximum . + StorableArray.mapShape Shape.Deferred) $ + mapAccumL + (\alphai xi -> + swap $ mapSnd (Vector.mul (emission hmm xi)) $ + matrixMaxMul (transition hmm) alphai) + (Vector.mul (emission hmm x) (initial hmm)) xs + +matrixMaxMul :: + (Shape.Indexed sh, Eq sh, Shape.Index sh ~ ix, Class.Real a) => + SquareMatrix sh a -> Vector sh a -> + (Vector (Shape.Deferred sh) (Shape.DeferredIndex ix), Vector sh a) +matrixMaxMul m v = + let sh = StorableArray.shape v + in mapPair (Vector.fromList (Shape.Deferred sh), Vector.fromList sh) $ + unzip $ + map (Vector.argAbsMaximum . + StorableArray.mapShape Shape.Deferred . + Vector.mul v) $ + Matrix.toRows m + + + +{- | +A trained model is a temporary form of a Hidden Markov model +that we need during the training on multiple training sequences. +It allows to collect knowledge over many sequences with 'mergeTrained', +even with mixed supervised and unsupervised training. +You finish the training by converting the trained model +back to a plain modul using 'finishTraining'. + +You can create a trained model in three ways: + +* supervised training using an emission sequence with associated states, + +* unsupervised training using an emission sequence and an existing Hidden Markov Model, + +* derive it from state sequence patterns, cf. "Math.HiddenMarkovModel.Pattern". +-} +data Trained distr sh prob = + Trained { + trainedInitial :: Vector sh prob, + trainedTransition :: SquareMatrix sh prob, + trainedDistribution :: distr + } + deriving (Show) + +instance + (NFData distr, NFData sh, NFData prob, Storable prob) => + NFData (Trained distr sh prob) where + rnf hmm = + rnf (trainedInitial hmm, trainedTransition hmm, trainedDistribution hmm) + + +sumTransitions :: + (Shape.C sh, Eq sh, Class.Real e) => + T distr sh e -> [SquareMatrix sh e] -> SquareMatrix sh e +sumTransitions hmm = + List.foldl' Vector.add + (Vector.constant (StorableArray.shape $ transition hmm) 0) + +{- | +Baum-Welch algorithm +-} +trainUnsupervised :: + (Distr.Estimate tdistr distr, Distr.StateShape distr ~ sh, Eq sh, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + T distr sh prob -> NonEmpty.T [] emission -> Trained tdistr sh prob +trainUnsupervised hmm xs = + let (recipLikelihood, alphas, betas) = alphaBeta hmm xs + zetas = zetaFromAlphaBeta recipLikelihood alphas betas + zeta0 = NonEmpty.head zetas + + in Trained { + trainedInitial = zeta0, + trainedTransition = + sumTransitions hmm $ + xiFromAlphaBeta hmm recipLikelihood xs alphas betas, + trainedDistribution = + Distr.accumulateEmissions $ + Array.fromList (StorableArray.shape zeta0) $ + map (zip (NonEmpty.flatten xs)) $ + List.transpose $ map Vector.toList $ NonEmpty.flatten zetas + } + + +mergeTrained :: + (Shape.C sh, Eq sh, + Distr.Estimate tdistr distr, Distr.Probability distr ~ prob) => + Trained tdistr sh prob -> Trained tdistr sh prob -> Trained tdistr sh prob +mergeTrained hmm0 hmm1 = + Trained { + trainedInitial = Vector.add (trainedInitial hmm0) (trainedInitial hmm1), + trainedTransition = + Vector.add (trainedTransition hmm0) (trainedTransition hmm1), + trainedDistribution = + Distr.combine + (trainedDistribution hmm0) (trainedDistribution hmm1) + } + +instance + (Shape.C sh, Eq sh, + Distr.Estimate tdistr distr, Distr.Probability distr ~ prob) => + Sg.Semigroup (Trained tdistr sh prob) where + (<>) = mergeTrained + + +toCells :: + (Distr.ToCSV distr, Shape.Indexed sh, Class.Real prob, Show prob) => + T distr sh prob -> [[String]] +toCells hmm = + (HMMCSV.cellsFromVector $ initial hmm) : + (HMMCSV.cellsFromSquare $ transition hmm) ++ + [] : + (Distr.toCells $ distribution hmm) + +parseCSV :: + (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh, Shape.C stateSh, + Class.Real prob, Read prob) => + (Int -> stateSh) -> HMMCSV.CSVParser (T distr stateSh prob) +parseCSV makeShape = do + v <- + StorableArray.mapShape (makeShape . Shape.zeroBasedSize) <$> + HMMCSV.parseNonEmptyVectorCells + let sh = StorableArray.shape v + m <- HMMCSV.parseSquareMatrixCells sh + HMMCSV.skipEmptyRow + distr <- Distr.parseCells sh + return $ Cons { + initial = v, + transition = m, + distribution = distr + } diff --git a/src/Math/HiddenMarkovModel/Test.hs b/src/Math/HiddenMarkovModel/Test.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Test.hs @@ -0,0 +1,258 @@ +{- | +Do not import this module, it is only intended for testing! +-} +module Math.HiddenMarkovModel.Test (tests) where + +import qualified Math.HiddenMarkovModel.Example.TrafficLightPrivate + as TrafficLight +import qualified Math.HiddenMarkovModel.Example.CirclePrivate as Circle + +import qualified Math.HiddenMarkovModel as HMM +import qualified Math.HiddenMarkovModel.Normalized as Normalized +import qualified Math.HiddenMarkovModel.Private as Priv +import qualified Math.HiddenMarkovModel.Distribution as Distr +import Math.HiddenMarkovModel.Utility (SquareMatrix, squareFromLists) + +import qualified Numeric.LAPACK.Vector as Vector +import qualified Numeric.LAPACK.ShapeStatic as ShapeStatic +import Numeric.LAPACK.Vector (Vector) + +import qualified Data.Array.Comfort.Shape as Shape + +import qualified Data.FixedLength as FL +import Data.FixedLength ((!:)) + +import qualified Type.Data.Num.Unary.Literal as TypeNum +import Type.Base.Proxy (Proxy(Proxy)) + +import qualified Test.QuickCheck as QC +import qualified System.Random as Rnd + +import Control.DeepSeq (deepseq) + +import qualified Data.NonEmpty.Class as NonEmptyC +import qualified Data.NonEmpty as NonEmpty +import qualified Data.Traversable as Trav +import qualified Data.Foldable as Fold +import qualified Data.Map as Map +import Data.Tuple.HT (mapSnd) + +import Text.Printf (printf) + + +type StateSet = ShapeStatic.ZeroBased TypeNum.U4 + +hmm :: HMM.Discrete Char StateSet Double +hmm = + HMM.Cons { + HMM.initial = stateVector 0.1 0.2 0.3 0.4, + HMM.transition = + squareFromLists stateSet $ + stateVector 0.7 0.1 0.0 0.2 : + stateVector 0.1 0.6 0.1 0.0 : + stateVector 0.1 0.2 0.7 0.0 : + stateVector 0.1 0.1 0.2 0.8 : + [], + HMM.distribution = + Distr.Discrete $ Map.fromList $ + ('a', stateVector 1 0 0 0) : + ('b', stateVector 0 1 0 1) : + ('c', stateVector 0 0 1 0) : + [] + } + +stateSet :: StateSet +stateSet = ShapeStatic.ZeroBased Proxy + +stateVector :: Double -> Double -> Double -> Double -> Vector StateSet Double +stateVector x0 x1 x2 x3 = ShapeStatic.vector $ x0!:x1!:x2!:x3!:FL.end +-- stateVector = FL.curry ShapeStatic.vector + + +sequ :: NonEmpty.T [] Char +sequ = NonEmpty.cons 'a' $ take 20 (HMM.generate hmm (Rnd.mkStdGen 42)) + +possibleStates :: Char -> [FL.Index TypeNum.U4] +possibleStates c = + map fst $ filter snd $ + zip (Shape.indices stateSet) $ + map + (\p -> + case p of + 0 -> False + 1 -> True + _ -> error "invalid emission probability (must be 0 or 1)") $ + Vector.toList $ + Map.findWithDefault (error "invalid character") c $ + case HMM.distribution hmm of Distr.Discrete m -> m + +{- | +Should all be equal. +-} +sequLikelihood :: ((Double, Double), Double, Double, NonEmpty.T [] Double) +sequLikelihood = + ((Priv.forward hmm sequ, Priv.backward hmm sequ), + exp $ Normalized.logLikelihood hmm sequ, + sum $ + map (NonEmpty.product . HMM.probabilitySequence hmm) $ + Trav.mapM (\c -> map (flip (,) c) $ possibleStates c) sequ, + NonEmptyC.zipWith Vector.dot + (Priv.alpha hmm sequ) + (Priv.beta hmm $ NonEmpty.tail sequ)) + +{- | +Should all be one. +-} +sequLikelihoodNormalized :: NonEmpty.T [] Double +sequLikelihoodNormalized = + let (calphas,betas) = Normalized.alphaBeta hmm sequ + in NonEmptyC.zipWith Vector.dot (fmap snd calphas) betas + + +{- | +Lists should be equal, but the first list contains one less element. +-} +zetas :: + ([Vector StateSet Double], + NonEmpty.T [] (Vector StateSet Double), + NonEmpty.T [] (Vector StateSet Double)) +zetas = + let (recipLikelihood, alphas, betas) = Priv.alphaBeta hmm sequ + in (Priv.zetaFromXi $ + Priv.xiFromAlphaBeta hmm recipLikelihood sequ alphas betas, + Priv.zetaFromAlphaBeta recipLikelihood alphas betas, + uncurry Normalized.zetaFromAlphaBeta $ + Normalized.alphaBeta hmm sequ) + + +distance :: + (Shape.C sh, Eq sh) => + Vector.Vector sh Double -> Vector.Vector sh Double -> Double +distance x y = Vector.normInf (Vector.sub x y) + + +{- | +Quick test of zetas - result should be @(True, very small, very small)@. +-} +zetasDiff :: (Bool, Double, Double) +zetasDiff = + case zetas of + (z0,z1,z2) -> + (length z0 == length (NonEmpty.tail z1) && + length z0 == length (NonEmpty.tail z2), + maximum $ zipWith distance z0 $ NonEmpty.init z1, + NonEmpty.maximum $ NonEmptyC.zipWith distance z1 z2) + +{- | +Lists should be equal +-} +xis :: ([SquareMatrix StateSet Double], [SquareMatrix StateSet Double]) +xis = + let (recipLikelihood, alphas, betas) = Priv.alphaBeta hmm sequ + in (Priv.xiFromAlphaBeta hmm recipLikelihood sequ alphas betas, + uncurry (Normalized.xiFromAlphaBeta hmm sequ) $ + Normalized.alphaBeta hmm sequ) + +{- | +Quick test of xis - result should be @(True, very small)@. +-} +xisDiff :: (Bool, Double) +xisDiff = + case xis of + (x0,x1) -> (length x0 == length x1, maximum $ zipWith distance x0 x1) + + +reveal :: Bool +reveal = + Normalized.reveal hmm sequ == Priv.reveal hmm sequ + + +trainUnsupervised :: + (HMM.DiscreteTrained Char StateSet Double, + HMM.DiscreteTrained Char StateSet Double) +trainUnsupervised = + (Priv.trainUnsupervised hmm sequ, + Normalized.trainUnsupervised hmm sequ) + +trainUnsupervisedDiff :: (Double, Double, (Bool, Double)) +trainUnsupervisedDiff = + case trainUnsupervised of + (hmm0,hmm1) -> + (distance (Priv.trainedTransition hmm0) (Priv.trainedTransition hmm1), + distance + (Priv.trainedInitial hmm0) (Priv.trainedInitial hmm1), + case (Priv.trainedDistribution hmm0, Priv.trainedDistribution hmm1) of + (Distr.DiscreteTrained m0, Distr.DiscreteTrained m1) -> + (Map.size m0 == Map.size m1, + Fold.maximum $ Map.intersectionWith distance m0 m1)) + + +nonEmptyScanr :: Int -> [Int] -> Bool +nonEmptyScanr x xs = + Normalized.nonEmptyScanr (-) x xs == NonEmpty.scanr (-) x xs + + +circleTraining :: (Int, Circle.HMM) -> Bool +circleTraining (maxDiff,hmm_) = + maxDiff >= + (length $ filter id $ NonEmpty.flatten $ + NonEmpty.zipWith (/=) + (HMM.reveal hmm_ Circle.circle) (fmap fst Circle.circleLabeled)) + + +allPair :: (a -> Bool, b -> Bool) -> (a,b) -> Bool +allPair (f,g) (a,b) = f a && g b + +allTriple :: (a -> Bool, b -> Bool, c -> Bool) -> (a,b,c) -> Bool +allTriple (f,g,h) (a,b,c) = f a && g b && h c + +almostZero :: Double -> Bool +almostZero x = x < 1e-10 + +almostOne :: Double -> Bool +almostOne x = almostZero $ abs (x-1) + +almostEqual :: Double -> Double -> Bool +almostEqual x y = almostZero $ abs (x-y) + +tests :: [(String, QC.Property)] +tests = + ("sequLikelihood", + QC.property $ + case sequLikelihood of + (forwardBackward, expLog, sumProb, alphaBetas) -> + allPair (almostEqual sumProb, almostEqual sumProb) forwardBackward + && + almostEqual sumProb expLog + && + length (NonEmpty.tail sequ) == length (NonEmpty.tail alphaBetas) + && + Fold.all (almostEqual sumProb) alphaBetas) : + ("sequLikelihoodNormalized", + QC.property $ + length (NonEmpty.tail sequ) == + length (NonEmpty.tail sequLikelihoodNormalized) + && + Fold.all almostOne sequLikelihoodNormalized) : + ("zetasDiff", + QC.property $ allTriple (id, almostZero, almostZero) zetasDiff) : + ("xisDiff", QC.property $ allPair (id, almostZero) xisDiff) : + ("reveal", QC.property reveal) : + ("trainUnsupervisedDiff", + QC.property $ + allTriple (almostZero, almostZero, allPair (id, almostZero)) $ + trainUnsupervisedDiff) : + ("nonEmptyScanr", QC.property nonEmptyScanr) : + (zip + (map (printf "TrafficLight.verifyRevelation.%d") [(0::Int) ..]) + (map QC.property TrafficLight.verifyRevelations)) ++ + ("TrafficLight.hmmIterativelyTrained.defined", + QC.property $ deepseq TrafficLight.hmmIterativelyTrained True) : + (map (mapSnd (QC.property . circleTraining)) $ + ("Circle.hmm", (0, Circle.hmm)) : + ("Circle.reconstructModel", (0, Circle.reconstructModel)) : + ("Circle.hmmTrainedSupervised", (0, Circle.hmmTrainedSupervised)) : + ("Circle.hmmTrainedUnsupervised", (0, Circle.hmmTrainedUnsupervised)) : + ("Circle.hmmIterativelyTrained", (40, Circle.hmmIterativelyTrained)) : + []) ++ + [] diff --git a/src/Math/HiddenMarkovModel/Utility.hs b/src/Math/HiddenMarkovModel/Utility.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Utility.hs @@ -0,0 +1,72 @@ +module Math.HiddenMarkovModel.Utility where + +import qualified Numeric.LAPACK.Matrix.Triangular as Triangular +import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian +import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape +import qualified Numeric.LAPACK.Matrix.Square as Square +import qualified Numeric.LAPACK.Matrix as Matrix +import qualified Numeric.LAPACK.Vector as Vector +import Numeric.LAPACK.Matrix.Triangular (Diagonal) +import Numeric.LAPACK.Vector (Vector) + +import qualified Numeric.Netlib.Class as Class + +import qualified Data.Array.Comfort.Storable as StorableArray +import qualified Data.Array.Comfort.Boxed as Array +import qualified Data.Array.Comfort.Shape as Shape + +import Foreign.Storable (Storable) + +import qualified System.Random as Rnd + +import qualified Control.Monad.Trans.State as MS + + +type SquareMatrix sh = Square.Square sh + +normalizeProb :: (Shape.C sh, Class.Real a) => Vector sh a -> Vector sh a +normalizeProb = snd . normalizeFactor + +normalizeFactor :: (Shape.C sh, Class.Real a) => Vector sh a -> (a, Vector sh a) +normalizeFactor xs = + let c = Vector.sum xs + in (c, Vector.scale (recip c) xs) + +-- see htam:Stochastic +randomItemProp :: + (Rnd.RandomGen g, Rnd.Random b, Num b, Ord b) => + [(a,b)] -> MS.State g a +randomItemProp props = + let (keys,ps) = unzip props + in do p <- MS.state (Rnd.randomR (0, sum ps)) + return $ + fst $ head $ dropWhile ((0<=) . snd) $ + zip keys $ tail $ scanl (-) p ps + +attachOnes :: (Num b) => [a] -> [(a,b)] +attachOnes = map (flip (,) 1) + + +vectorDim :: Shape.C sh => Vector sh a -> Int +vectorDim = Shape.size . StorableArray.shape + +singleton :: (Class.Real a) => a -> Vector () a +singleton = Vector.constant () + + +hermitianFromList :: + (Shape.C sh, Storable a) => sh -> [a] -> Hermitian.Hermitian sh a +hermitianFromList = Hermitian.fromList MatrixShape.RowMajor + + +squareConstant :: + (Shape.C sh, Class.Real a) => sh -> a -> SquareMatrix sh a +squareConstant = Vector.constant . MatrixShape.square MatrixShape.RowMajor + +squareFromLists :: + (Shape.C sh, Eq sh, Storable a) => sh -> [Vector sh a] -> SquareMatrix sh a +squareFromLists sh = + Square.fromGeneral . Matrix.fromRowArray sh . Array.fromList sh + +diagonal :: (Shape.C sh, Class.Real a) => Vector sh a -> Diagonal sh a +diagonal = Triangular.diagonal MatrixShape.RowMajor diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Main.hs @@ -0,0 +1,10 @@ +module Main where + +import Math.HiddenMarkovModel.Test (tests) + +import qualified Test.QuickCheck as QC + + +main :: IO () +main = + mapM_ (\(name,prop) -> putStr (name ++ ": ") >> QC.quickCheck prop) tests