diff --git a/Changes.md b/Changes.md --- a/Changes.md +++ b/Changes.md @@ -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. diff --git a/hmm-lapack.cabal b/hmm-lapack.cabal --- a/hmm-lapack.cabal +++ b/hmm-lapack.cabal @@ -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 diff --git a/private/Math/HiddenMarkovModel/CSV.hs b/private/Math/HiddenMarkovModel/CSV.hs new file mode 100644 --- /dev/null +++ b/private/Math/HiddenMarkovModel/CSV.hs @@ -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) diff --git a/private/Math/HiddenMarkovModel/Example/CirclePrivate.hs b/private/Math/HiddenMarkovModel/Example/CirclePrivate.hs new file mode 100644 --- /dev/null +++ b/private/Math/HiddenMarkovModel/Example/CirclePrivate.hs @@ -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 diff --git a/private/Math/HiddenMarkovModel/Example/SineWavePrivate.hs b/private/Math/HiddenMarkovModel/Example/SineWavePrivate.hs new file mode 100644 --- /dev/null +++ b/private/Math/HiddenMarkovModel/Example/SineWavePrivate.hs @@ -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 diff --git a/private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs b/private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs new file mode 100644 --- /dev/null +++ b/private/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs @@ -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) diff --git a/private/Math/HiddenMarkovModel/Normalized.hs b/private/Math/HiddenMarkovModel/Normalized.hs new file mode 100644 --- /dev/null +++ b/private/Math/HiddenMarkovModel/Normalized.hs @@ -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 + } diff --git a/private/Math/HiddenMarkovModel/Private.hs b/private/Math/HiddenMarkovModel/Private.hs new file mode 100644 --- /dev/null +++ b/private/Math/HiddenMarkovModel/Private.hs @@ -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 + } diff --git a/private/Math/HiddenMarkovModel/Public.hs b/private/Math/HiddenMarkovModel/Public.hs new file mode 100644 --- /dev/null +++ b/private/Math/HiddenMarkovModel/Public.hs @@ -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 diff --git a/private/Math/HiddenMarkovModel/Public/Distribution.hs b/private/Math/HiddenMarkovModel/Public/Distribution.hs new file mode 100644 --- /dev/null +++ b/private/Math/HiddenMarkovModel/Public/Distribution.hs @@ -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)) diff --git a/private/Math/HiddenMarkovModel/Utility.hs b/private/Math/HiddenMarkovModel/Utility.hs new file mode 100644 --- /dev/null +++ b/private/Math/HiddenMarkovModel/Utility.hs @@ -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) diff --git a/src/Math/HiddenMarkovModel.hs b/src/Math/HiddenMarkovModel.hs --- a/src/Math/HiddenMarkovModel.hs +++ b/src/Math/HiddenMarkovModel.hs @@ -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 diff --git a/src/Math/HiddenMarkovModel/CSV.hs b/src/Math/HiddenMarkovModel/CSV.hs deleted file mode 100644 --- a/src/Math/HiddenMarkovModel/CSV.hs +++ /dev/null @@ -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) diff --git a/src/Math/HiddenMarkovModel/Distribution.hs b/src/Math/HiddenMarkovModel/Distribution.hs --- a/src/Math/HiddenMarkovModel/Distribution.hs +++ b/src/Math/HiddenMarkovModel/Distribution.hs @@ -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 diff --git a/src/Math/HiddenMarkovModel/Example/CirclePrivate.hs b/src/Math/HiddenMarkovModel/Example/CirclePrivate.hs deleted file mode 100644 --- a/src/Math/HiddenMarkovModel/Example/CirclePrivate.hs +++ /dev/null @@ -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 diff --git a/src/Math/HiddenMarkovModel/Example/SineWave.hs b/src/Math/HiddenMarkovModel/Example/SineWave.hs --- a/src/Math/HiddenMarkovModel/Example/SineWave.hs +++ b/src/Math/HiddenMarkovModel/Example/SineWave.hs @@ -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 diff --git a/src/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs b/src/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs deleted file mode 100644 --- a/src/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs +++ /dev/null @@ -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) diff --git a/src/Math/HiddenMarkovModel/Normalized.hs b/src/Math/HiddenMarkovModel/Normalized.hs deleted file mode 100644 --- a/src/Math/HiddenMarkovModel/Normalized.hs +++ /dev/null @@ -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 - } diff --git a/src/Math/HiddenMarkovModel/Pattern.hs b/src/Math/HiddenMarkovModel/Pattern.hs --- a/src/Math/HiddenMarkovModel/Pattern.hs +++ b/src/Math/HiddenMarkovModel/Pattern.hs @@ -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 diff --git a/src/Math/HiddenMarkovModel/Private.hs b/src/Math/HiddenMarkovModel/Private.hs deleted file mode 100644 --- a/src/Math/HiddenMarkovModel/Private.hs +++ /dev/null @@ -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 - } diff --git a/src/Math/HiddenMarkovModel/Test.hs b/src/Math/HiddenMarkovModel/Test.hs deleted file mode 100644 --- a/src/Math/HiddenMarkovModel/Test.hs +++ /dev/null @@ -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)) : - []) ++ - [] diff --git a/src/Math/HiddenMarkovModel/Utility.hs b/src/Math/HiddenMarkovModel/Utility.hs deleted file mode 100644 --- a/src/Math/HiddenMarkovModel/Utility.hs +++ /dev/null @@ -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) diff --git a/test/Main.hs b/test/Main.hs --- a/test/Main.hs +++ b/test/Main.hs @@ -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 diff --git a/test/Test.hs b/test/Test.hs new file mode 100644 --- /dev/null +++ b/test/Test.hs @@ -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) : + [] diff --git a/test/Test/Main.hs b/test/Test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Test/Main.hs @@ -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 diff --git a/test/Test/Math/HiddenMarkovModel/Example/CirclePrivate.hs b/test/Test/Math/HiddenMarkovModel/Example/CirclePrivate.hs new file mode 100644 --- /dev/null +++ b/test/Test/Math/HiddenMarkovModel/Example/CirclePrivate.hs @@ -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)) diff --git a/test/Test/Math/HiddenMarkovModel/Example/SineWavePrivate.hs b/test/Test/Math/HiddenMarkovModel/Example/SineWavePrivate.hs new file mode 100644 --- /dev/null +++ b/test/Test/Math/HiddenMarkovModel/Example/SineWavePrivate.hs @@ -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]"]] diff --git a/test/Test/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs b/test/Test/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs new file mode 100644 --- /dev/null +++ b/test/Test/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs @@ -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) diff --git a/test/Test/Math/HiddenMarkovModel/Normalized.hs b/test/Test/Math/HiddenMarkovModel/Normalized.hs new file mode 100644 --- /dev/null +++ b/test/Test/Math/HiddenMarkovModel/Normalized.hs @@ -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]))