hmm-lapack 0.3.0.3 → 0.4

raw patch · 13 files changed

+595/−564 lines, 13 filesdep −boxesdep ~comfort-arraydep ~lapackdep ~non-emptyPVP ok

version bump matches the API change (PVP)

Dependencies removed: boxes

Dependency ranges changed: comfort-array, lapack, non-empty, random, semigroups

API changes (from Hackage documentation)

Files

• hmm-lapack.cabal +6/−7
• src/Math/HiddenMarkovModel.hs +43/−58
• src/Math/HiddenMarkovModel/CSV.hs +9/−9
• src/Math/HiddenMarkovModel/Distribution.hs +320/−258
• src/Math/HiddenMarkovModel/Example/CirclePrivate.hs +1/−1
• src/Math/HiddenMarkovModel/Example/SineWave.hs +1/−1
• src/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs +5/−11
• src/Math/HiddenMarkovModel/Named.hs +15/−18
• src/Math/HiddenMarkovModel/Normalized.hs +33/−51
• src/Math/HiddenMarkovModel/Pattern.hs +24/−22
• src/Math/HiddenMarkovModel/Private.hs +107/−111
• src/Math/HiddenMarkovModel/Test.hs +17/−10
• src/Math/HiddenMarkovModel/Utility.hs +14/−7

@@ -1,5 +1,5 @@ Name:                hmm-lapack-Version:             0.3.0.3+Version:             0.4 Synopsis:            Hidden Markov Models using LAPACK primitives Description:   Hidden Markov Models implemented using LAPACK data types and operations.@@ -41,7 +41,7 @@   Changes.md  Source-Repository this-  Tag:         0.3.0.3+  Tag:         0.4   Type:        darcs   Location:    http://hub.darcs.net/thielema/hmm-lapack @@ -67,19 +67,18 @@     Math.HiddenMarkovModel.Utility     Math.HiddenMarkovModel.CSV   Build-Depends:-    lapack >=0.2.2 && <0.3,+    lapack >=0.3 && <0.4,     fixed-length >=0.2.1 && <0.3,     tfp >=1.0 && <1.1,     netlib-ffi >=0.1.1 && <0.2,-    comfort-array >=0.2 && <0.4,+    comfort-array >=0.4 && <0.5,     QuickCheck >=2.5 && <3,     explicit-exception >=0.1.7 && <0.2,-    boxes >=0.1.5 && <0.2,     lazy-csv >=0.5 && <0.6,     random >=1.0 && <1.2,     transformers >= 0.2 && <0.6,-    non-empty >=0.2.1 && <0.4,-    semigroups >=0.17 && <0.19,+    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,

@@ -26,17 +26,18 @@ import Math.HiddenMarkovModel.Private           (T(..), Trained(..), mergeTrained, toCells, parseCSV) import Math.HiddenMarkovModel.Utility-          (SquareMatrix, squareConstant, distance,+          (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 Data.Array.Comfort.Boxed as Array  import qualified Text.CSV.Lazy.String as CSV @@ -53,12 +54,12 @@   type DiscreteTrained symbol sh prob =-         Trained (Distr.DiscreteTrained symbol sh prob) sh prob-type Discrete symbol sh prob = T (Distr.Discrete symbol sh prob) sh prob+         Trained (Distr.Discrete symbol) sh prob+type Discrete symbol sh prob = T (Distr.Discrete symbol) sh prob  type GaussianTrained emiSh stateSh a =-         Trained (Distr.GaussianTrained emiSh stateSh a) stateSh a-type Gaussian emiSh stateSh a = T (Distr.Gaussian emiSh stateSh a) stateSh a+         Trained (Distr.Gaussian emiSh) stateSh a+type Gaussian emiSh stateSh a = T (Distr.Gaussian emiSh) stateSh a   {- |@@ -68,9 +69,8 @@ You can use this as a starting point for 'Normalized.trainUnsupervised'. -} uniform ::-   (Distr.Info distr, Distr.StateShape distr ~ sh, Shape.C sh,-    Distr.Probability distr ~ prob) =>-   distr -> T distr sh prob+   (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@@ -82,31 +82,29 @@   probabilitySequence ::-   (Traversable f, Distr.EmissionProb distr,-    Distr.StateShape distr ~ sh, Shape.Indexed sh, Shape.Index sh ~ state,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>-   T distr sh prob -> f (state, emission) -> f prob+   (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 StorableArray.!) . flip (,) s,+         ((transition hmm #!) . flip (,) s,           index s * Distr.emissionStateProb (distribution hmm) e s))       (initial hmm StorableArray.!)  generate ::-   (Rnd.RandomGen g, Ord prob, Rnd.Random prob, Distr.Generate distr,-    Distr.StateShape distr ~ sh, Shape.Indexed sh, Shape.Index sh ~ state,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>-   T distr sh prob -> g -> [emission]+   (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 ::-   (Rnd.RandomGen g, Ord prob, Rnd.Random prob, Distr.Generate distr,-    Distr.StateShape distr ~ sh, Shape.Indexed sh, Shape.Index sh ~ state,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>-   T distr sh prob -> g -> [(state, emission)]+   (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) $@@ -123,31 +121,24 @@ Contribute a manually labeled emission sequence to a HMM training. -} trainSupervised ::-   (Distr.StateShape distr ~ sh, Shape.Index sh ~ state,-    Distr.Estimate tdistr distr,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>-   sh -> NonEmpty.T [] (state, emission) -> Trained tdistr sh prob+   (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 =-             StorableArray.fromAssociations sh 0-                [(getState (NonEmpty.head xs), 1)],+          trainedInitial = Vector.unit sh $ getState $ NonEmpty.head xs,           trainedTransition =-             Matrix.transpose $-             StorableArray.accumulate (+) (squareConstant sh 0) $+             Matrix.transpose $ ArrMatrix.fromVector $+             StorableArray.accumulate (+)+                (ArrMatrix.toVector $ squareConstant sh 0) $              attachOnes $ NonEmpty.mapAdjacent (,) $ fmap getState xs,-          trainedDistribution =-             Distr.accumulateEmissions $ Array.map attachOnes $-             Array.accumulate (flip (:))-                (Array.fromList sh $ replicate (Shape.size sh) [])-                (NonEmpty.flatten xs)+          trainedDistribution = Distr.accumulateEmissions sh xs        }  finishTraining ::-   (Shape.C sh, Eq sh,-    Distr.Estimate tdistr distr, Distr.Probability distr ~ prob) =>-   Trained tdistr sh prob -> T distr sh prob+   (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,@@ -156,18 +147,15 @@    }  normalizeProbColumns ::-   (Shape.C sh, Eq sh, Class.Real a) => SquareMatrix sh a -> SquareMatrix sh a+   (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 ::-   (Shape.C sh, Eq sh,-    Distr.Estimate tdistr distr, Distr.Probability distr ~ prob,-    Foldable f) =>-   (trainingData -> Trained tdistr sh prob) ->-   NonEmpty.T f trainingData -> T distr sh prob-trainMany train =-   finishTraining . NonEmpty.foldl1Map mergeTrained train+   (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   @@ -181,25 +169,22 @@ should suffice for defining an abort criterion. -} deviation ::-   (Shape.C sh, Eq sh, Class.Real prob, Ord prob) =>-   T distr sh prob -> T distr sh prob -> prob+   (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`-   distance (transition hmm0) (transition hmm1)+   matrixDistance (transition hmm0) (transition hmm1)   toCSV ::-   (Distr.ToCSV distr, Shape.Indexed sh, Class.Real prob, Show prob) =>-   T distr sh prob -> String+   (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+   CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $ toCells hmm  fromCSV ::-   (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh,-    Shape.Indexed stateSh, Shape.Index stateSh ~ state,-    Class.Real prob, Read prob) =>-   (Int -> stateSh) -> String -> ME.Exceptional String (T distr stateSh prob)+   (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

@@ -1,16 +1,15 @@ module Math.HiddenMarkovModel.CSV where -import Math.HiddenMarkovModel.Utility (SquareMatrix, vectorDim)+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 (ZeroInt)+import Numeric.LAPACK.Matrix (ShapeInt) import Numeric.LAPACK.Vector (Vector)  import qualified Numeric.Netlib.Class as Class -import qualified Data.Array.Comfort.Storable as ComfortArray import qualified Data.Array.Comfort.Shape as Shape  import qualified Text.CSV.Lazy.String as CSV@@ -32,7 +31,7 @@ cellsFromVector = map show . Vector.toList  cellsFromSquare ::-   (Shape.Indexed sh, Show a, Class.Real a) => SquareMatrix sh a -> [[String]]+   (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]]@@ -110,20 +109,21 @@  parseVectorCells ::    (Read a, Class.Real a) =>-   CSVParser (Vector ZeroInt 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 ZeroInt 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 ZeroInt a)+   CSVParser (Vector ShapeInt a) parseNonEmptyVectorCells = do    v <- parseVectorCells    assert (vectorDim v > 0) "no data for vector"@@ -143,14 +143,14 @@  parseSquareMatrixCells ::    (Shape.C sh, Read a, Class.Real a) =>-   sh -> CSVParser (SquareMatrix sh 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 $-      ComfortArray.reshape (MatrixShape.square MatrixShape.RowMajor sh) $+      Matrix.reshape (MatrixShape.square MatrixShape.RowMajor sh) $       Matrix.fromRows (Shape.ZeroBased n) rows  parseStringList :: CSV.CSVRow -> CSVParser [String]

@@ -1,28 +1,34 @@ {-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE EmptyDataDecls #-} module Math.HiddenMarkovModel.Distribution (-   Emission, Probability, StateShape,-   Info(..), Generate(..), EmissionProb(..), Estimate(..),+   T(..), Trained(..), Emission,+   Show(..), NFData(..), Format(..),+   Info(..), Generate(..), EmissionProb(..),+   Estimate(..), accumulateEmissionVectors, -   Discrete(..), DiscreteTrained(..),-   Gaussian(..), GaussianTrained(..), gaussian,+   Discrete, discreteFromList,+   Gaussian, gaussian, gaussianTrained,     ToCSV(..), FromCSV(..), HMMCSV.CSVParser, CSVSymbol(..),    ) where  import qualified Math.HiddenMarkovModel.CSV as HMMCSV-import Math.HiddenMarkovModel.Utility (SquareMatrix, randomItemProp, vectorDim)+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 Numeric.LAPACK.Matrix ((<#))+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, Format(format))+import Numeric.LAPACK.Format (FormatArray)+import Numeric.LAPACK.Output (Output)  import qualified Numeric.Netlib.Class as Class import Foreign.Storable (Storable)@@ -31,334 +37,397 @@ 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 qualified Text.PrettyPrint.Boxes as TextBox-import Text.PrettyPrint.Boxes ((<>), (<+>)) import Text.Read.HT (maybeRead) import Text.Printf (printf)  import qualified Control.Monad.Exception.Synchronous as ME import qualified Control.Monad.Trans.Class as MT import qualified Control.Monad.Trans.State as MS-import Control.DeepSeq (NFData, rnf)+import qualified Control.DeepSeq as DeepSeq import Control.Monad (liftM2)-import Control.Applicative (liftA2, (<|>))+import Control.Applicative (liftA2) +import qualified Data.NonEmpty.Map as NonEmptyMap import qualified Data.NonEmpty as NonEmpty-import qualified Data.Foldable as Fold+import qualified Data.Semigroup as Sg import qualified Data.Map as Map-import qualified Data.Set as Set import qualified Data.List.HT as ListHT import qualified Data.List as List import Data.Functor.Identity (Identity(Identity), runIdentity)-import Data.Foldable (Foldable, foldMap)-import Data.Tuple.HT (mapFst, fst3, swap)-import Data.Monoid (Endo(Endo, appEndo))-import Data.Map (Map)-import Data.Maybe (fromMaybe, listToMaybe)+import Data.Tuple.HT (snd3)+import Data.Set (Set)+import Data.Maybe (listToMaybe) -import Prelude ()-import Prelude2010+import qualified Prelude as P+import Prelude2010 hiding (Show, showsPrec)  -type HermitianMatrix sh = Hermitian.Hermitian sh-type UpperTriangular sh = Triangular.Upper sh +data family T typ sh prob+data family Trained typ sh prob -type family Probability distr-type family Emission distr-type family StateShape distr+type family Emission typ prob  -class (Class.Real (Probability distr)) => Info distr where-   statesShape :: distr -> StateShape distr+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 -class (Class.Real (Probability distr)) => Generate distr where+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 ::-      (Rnd.RandomGen g, Emission distr ~ emission, StateShape distr ~ sh) =>-      distr -> Shape.Index sh -> MS.State g emission+      (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-   (Shape.Indexed (StateShape distr), Class.Real (Probability distr)) =>-      EmissionProb distr where+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 ::-      distr -> Emission distr -> Vector (StateShape distr) (Probability distr)+      (Shape.C sh, Class.Real prob) =>+      T typ sh prob -> Emission typ prob -> Vector sh prob    emissionStateProb ::-      distr -> Emission distr -> Shape.Index (StateShape distr) -> Probability distr+      (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-   (Distribution tdistr ~ distr, Trained distr ~ tdistr, EmissionProb distr) =>-      Estimate tdistr distr where-   type Distribution tdistr-   type Trained distr+class (EmissionProb typ) => Estimate typ where    accumulateEmissions ::-      (Probability distr ~ prob, StateShape distr ~ sh) =>-      Array sh [(Emission distr, prob)] -> tdistr-   -- could as well be in Semigroup class-   combine :: tdistr -> tdistr -> tdistr-   normalize :: tdistr -> distr+      (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) -newtype Discrete symbol sh prob = Discrete (Map symbol (Vector sh prob))-   deriving (Show)+instance+   (Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>+      Sg.Semigroup (Trained typ sh prob) where+   (<>) = combine -newtype-   DiscreteTrained symbol sh prob =-      DiscreteTrained (Map symbol (Vector sh prob))-   deriving (Show) -type instance Probability (Discrete symbol sh prob) = prob-type instance Emission (Discrete symbol sh prob) = symbol-type instance StateShape (Discrete symbol sh prob) = sh+data Discrete symbol +newtype instance T (Discrete symbol) sh prob =+      Discrete (Matrix.General (Set symbol) sh prob) -instance-   (NFData sh, NFData prob, NFData symbol) =>-      NFData (Discrete symbol sh prob) where-   rnf (Discrete m) = rnf m+newtype instance Trained (Discrete symbol) sh prob =+      DiscreteTrained (NonEmptyMap.T symbol (Vector sh prob)) -instance-   (NFData sh, NFData prob, NFData symbol) =>-      NFData (DiscreteTrained symbol sh prob) where-   rnf (DiscreteTrained m) = rnf m+type instance Emission (Discrete symbol) prob = symbol -instance-   (FormatArray sh, Class.Real prob, Format symbol) =>-      Format (Discrete symbol sh prob) where++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) =-      TextBox.vsep 1 TextBox.left $-      map (\(sym,v) -> format fmt sym <> TextBox.char ':' <+> format fmt v) $-      Map.toAscList m+      Output.formatAligned $+      map (\(sym,v) ->+            map (Identity . Output.text) $+            (show sym ++ ":") : map (printFmt fmt) (Vector.toList v)) $+      Array.toAssociations $ Matrix.toRowArray m -instance-   (Shape.C sh, Class.Real prob, Ord symbol) =>-      Info (Discrete symbol sh prob) where-   statesShape (Discrete m) = StorableArray.shape $ snd $ Map.findMin m+-- cf. Data.Bifunctor.Flip+newtype Flip f b a = Flip {getFlip :: f a b} -instance-   (Shape.Indexed sh, Class.Real prob, Ord symbol, Ord prob, Rnd.Random prob) =>-      Generate (Discrete symbol sh prob) where-   generate (Discrete m) state =-      randomItemProp $ Map.toAscList $ fmap (StorableArray.! state) m+printFmt :: (Class.Real a) => String -> a -> String+printFmt fmt =+   getFlip $ Class.switchReal (Flip $ printf fmt) (Flip $ printf fmt) -instance-   (Shape.Indexed sh, Class.Real prob, Ord symbol) =>-      EmissionProb (Discrete symbol sh prob) where-   emissionProb (Discrete m) =-      mapLookup "emitDiscrete: unknown emission symbol" m+instance (Ord symbol) => Info (Discrete symbol) where+   statesShape (Discrete m) = Matrix.width m+   statesShapeTrained (DiscreteTrained m) = discreteStateShape m -instance-   (Shape.Indexed sh, Eq sh, Class.Real prob, Ord symbol) =>-      Estimate (DiscreteTrained symbol sh prob) (Discrete symbol sh prob) where-   type Distribution (DiscreteTrained symbol sh prob) = Discrete symbol sh prob-   type Trained (Discrete symbol sh prob) = DiscreteTrained symbol sh prob-   accumulateEmissions grouped =-      let set = Set.toAscList $ foldMap (Set.fromList . map fst) grouped-          emi = Map.fromAscList $ zip set [0..]-      in  DiscreteTrained $ Map.fromAscList $ zip set $-          transposeVectorList $-          Array.map-             (StorableArray.accumulate (+)-                 (Vector.constant (Shape.ZeroBased $ length set) 0) .-              map (mapFst-                 (mapLookup "estimateDiscrete: unknown emission symbol" emi)))-             grouped+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 $ Map.unionWith Vector.add distr0 distr1+      DiscreteTrained $ NonEmptyMap.unionWith Vector.add distr0 distr1    normalize (DiscreteTrained distr) =-      Discrete $ if Map.null distr then distr else normalizeProbVecs distr+      Discrete $ normalizeProbColumns $ discreteFromMap distr -transposeVectorList ::-   (Shape.C sh, Eq sh, Class.Real a) =>-   Array sh (Vector Matrix.ZeroInt a) -> [Vector sh a]-transposeVectorList xs =-   case Array.toList xs of-      [] -> []-      x:_ -> Matrix.toRows $ Matrix.fromColumnArray (StorableArray.shape x) xs+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 -normalizeProbVecs ::-   (Shape.C sh, Eq sh, Foldable f, Functor f, Class.Real a) =>-   f (Vector sh a) -> f (Vector sh a)-normalizeProbVecs vs =-   let factors =-         StorableArray.map recip $ List.foldl1' Vector.add $ Fold.toList vs-   in fmap (Vector.mul factors) vs+discreteStateShape ::+   (Shape.C sh) => NonEmptyMap.T symbol (Vector sh prob) -> sh+discreteStateShape =+   StorableArray.shape . snd . fst . NonEmptyMap.minViewWithKey -mapLookup :: (Ord k) => String -> Map.Map k a -> k -> a-mapLookup msg dict x = Map.findWithDefault (error msg) x dict+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 -newtype Gaussian emiSh stateSh a =-      Gaussian (Array stateSh (Vector emiSh a, UpperTriangular emiSh a, a))-   deriving (Show) -newtype GaussianTrained emiSh stateSh a =++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-      (Array stateSh-         (Maybe (Vector emiSh a, HermitianMatrix emiSh a, a)))-   deriving (Show)+      (StorableArray.Array (stateSh, MatrixShape.Hermitian (():+:emiSh)) a) -type instance Probability (Gaussian emiSh stateSh a) = a-type instance Emission (Gaussian emiSh stateSh a) = Vector emiSh a-type instance StateShape (Gaussian emiSh stateSh a) = stateSh+type instance Emission (Gaussian emiSh) a = Vector emiSh a  -instance-   (NFData emiSh, NFData stateSh, Shape.C stateSh, NFData a, Storable a) =>-      NFData (Gaussian emiSh stateSh a) where-   rnf (Gaussian params) = rnf params+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-   (NFData emiSh, NFData stateSh, Shape.C stateSh, NFData a, Storable a) =>-      NFData (GaussianTrained emiSh stateSh a) where-   rnf (GaussianTrained params) = rnf params+instance (DeepSeq.NFData emiSh) => NFData (Gaussian emiSh) where+   rnf (Gaussian params) = DeepSeq.rnf params+   rnfTrained (GaussianTrained params) = DeepSeq.rnf params  -instance-   (FormatArray emiSh, Shape.C stateSh, Class.Real a) =>-      Format (Gaussian emiSh stateSh a) where+instance (FormatArray emiSh) => Format (Gaussian emiSh) where    format = runFormatGaussian $ Class.switchReal formatGaussian formatGaussian -newtype FormatGaussian emiSh stateSh a =+newtype FormatGaussian out emiSh stateSh a =    FormatGaussian-      {runFormatGaussian :: String -> Gaussian emiSh stateSh a -> TextBox.Box}+      {runFormatGaussian :: String -> T (Gaussian emiSh) stateSh a -> out}  formatGaussian ::-   (FormatArray emiSh, Shape.C stateSh, Class.Real a, Format a) =>-   FormatGaussian emiSh stateSh a+   (FormatArray emiSh, Shape.C stateSh,+    Class.Real a, Format.Format a, Output out) =>+   FormatGaussian out emiSh stateSh a formatGaussian =-   FormatGaussian $ \fmt (Gaussian params) -> format fmt $ Array.toList params+   FormatGaussian $ \fmt (Gaussian params) ->+      Format.format fmt $ Array.toList params  -instance-   (Shape.Indexed stateSh, Eq stateSh, Class.Real a) =>-      Info (Gaussian emiSh stateSh a) where+instance Info (Gaussian emiSh) where    statesShape (Gaussian params) = Array.shape params+   statesShapeTrained (GaussianTrained params) =+      fst $ StorableArray.shape params -instance-   (Shape.C emiSh, Eq emiSh, Shape.Indexed stateSh, Eq stateSh, Class.Real a) =>-      Generate (Gaussian emiSh stateSh a) where+instance (Shape.C emiSh, Eq emiSh) => Generate (Gaussian emiSh) where    generate (Gaussian allParams) state = do-      let (center, covarianceChol, _c) = allParams ! state+      let (_c, center, covarianceChol) = allParams ! state       seed <- MS.state Rnd.random       return $          Vector.add center $          Vector.random Vector.Normal (StorableArray.shape center) seed-            <# covarianceChol+            -*# covarianceChol -instance-   (Shape.C emiSh, Eq emiSh, Shape.Indexed stateSh, Eq stateSh, Class.Real a) =>-      EmissionProb (Gaussian emiSh stateSh a) where+instance (Shape.C emiSh, Eq emiSh) => EmissionProb (Gaussian emiSh) where+   mapStatesShape f (Gaussian m) = Gaussian $ Array.mapShape f m    emissionProb (Gaussian allParams) x =-      Vector.fromList (Array.shape allParams) $-      map (emissionProbGen x) $ Array.toList allParams+      StorableArray.fromBoxed $ fmap (gaussianEmissionProb x) allParams    emissionStateProb (Gaussian allParams) x s =-      emissionProbGen x $ allParams ! s+      gaussianEmissionProb x $ allParams ! s -emissionProbGen ::+gaussianEmissionProb ::    (Shape.C emiSh, Eq emiSh, Class.Real a) =>-   Vector emiSh a -> (Vector emiSh a, UpperTriangular emiSh a, a) -> a-emissionProbGen x (center, covarianceChol, c) =-   let x0 =-         Matrix.solveVector (Triangular.transpose covarianceChol) $-         Vector.sub x center-   in  c * cexp ((-1/2) * Vector.inner x0 x0)+   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) -instance-   (Shape.C emiSh, Eq emiSh, Shape.Indexed stateSh, Eq stateSh, Class.Real a) =>-      Estimate-         (GaussianTrained emiSh stateSh a)-         (Gaussian emiSh stateSh a) where-   type Distribution (GaussianTrained emiSh stateSh a) =-            Gaussian emiSh stateSh a-   type Trained (Gaussian emiSh stateSh a) = GaussianTrained emiSh stateSh a-   accumulateEmissions =-      let params xs =-            (NonEmpty.foldl1Map Vector.add (uncurry $ flip Vector.scale) xs,-             covarianceReal $ fmap swap xs,-             Fold.sum $ fmap snd xs)-      in  GaussianTrained . fmap (fmap params . NonEmpty.fetch)-   combine (GaussianTrained distr0) (GaussianTrained distr1) =-      let comb (center0, covariance0, weight0)-               (center1, covariance1, weight1) =-             (Vector.add center0 center1,-              Vector.add covariance0 covariance1,-              weight0 + weight1)-      in  GaussianTrained $ Array.zipWith (maybePlus comb) distr0 distr1+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 distr) =-      let params (centerSum, covarianceSum, weight) =-             let c = recip weight-                 center = Vector.scale c centerSum+   normalize (GaussianTrained m) =+      let params (weight, centerSum, covarianceSum) =+             let c = recip (weight#!((),()))+                 center = Vector.scale c $ Matrix.flattenRow centerSum              in  (center,-                  Vector.sub (Vector.scale c covarianceSum)+                  Matrix.sub+                     (Matrix.scaleRealReal c covarianceSum)                      (Hermitian.outer MatrixShape.RowMajor center))-      in  Gaussian $-          fmap-             (gaussianParameters . params .-              fromMaybe-                (error "Distribution.normalize: undefined array element")) $-          distr--maybePlus :: (a -> a -> a) -> Maybe a -> Maybe a -> Maybe a-maybePlus f mx my = liftA2 f mx my <|> mx <|> my+      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 ()) -newtype CovarianceReal f emiSh a =-   CovarianceReal-      {getCovarianceReal :: f (a, Vector emiSh a) -> HermitianMatrix emiSh a}+{- |+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) -covarianceReal ::-   (Shape.C emiSh, Eq emiSh, Class.Real a) =>-   NonEmpty.T [] (a, Vector emiSh a) -> HermitianMatrix emiSh a-covarianceReal =-   getCovarianceReal $-   Class.switchReal-      (CovarianceReal $ Hermitian.sumRank1NonEmpty MatrixShape.RowMajor)-      (CovarianceReal $ Hermitian.sumRank1NonEmpty MatrixShape.RowMajor)+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, HermitianMatrix emiSh prob) ->-   Gaussian emiSh stateSh prob-gaussian = consGaussian . fmap gaussianParameters+   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, HermitianMatrix emiSh prob) ->-   (Vector emiSh prob, UpperTriangular emiSh prob, prob)+   (Vector emiSh prob, Matrix.Hermitian emiSh prob) ->+   (prob, Vector emiSh prob, Triangular.Upper emiSh prob) gaussianParameters (center, covariance) =    gaussianFromCholesky center $ HermitianPD.decompose covariance -consGaussian ::-   (Shape.C stateSh) =>-   Array stateSh (Vector emiSh a, UpperTriangular emiSh a, a) ->-   Gaussian emiSh stateSh a-consGaussian = Gaussian- gaussianFromCholesky ::    (Shape.C emiSh, Class.Real prob) =>-   Vector emiSh prob -> UpperTriangular emiSh prob ->-   (Vector emiSh prob, UpperTriangular emiSh prob, prob)+   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  (center, covarianceChol,-        recip (sqrt2pi ^ vectorDim center * covarianceSqrtDet))+   in  (recip (sqrt2pi ^ vectorDim center * covarianceSqrtDet),+        center, covarianceChol)  sqrt2pi :: (Class.Real a) => a sqrt2pi = runIdentity $ Class.switchReal sqrt2piAux sqrt2piAux@@ -366,16 +435,16 @@ sqrt2piAux :: (Floating a) => Identity a sqrt2piAux = Identity $ sqrt (2*pi) -cexp :: (Class.Real a) => a -> a-cexp = appEndo $ Class.switchReal (Endo exp) (Endo exp) ---class ToCSV distr where-   toCells :: distr -> [[String]]+class ToCSV typ where+   toCells ::+      (Shape.C sh, Class.Real prob, P.Show prob) =>+      T typ sh prob -> [[String]] -class FromCSV distr where-   parseCells :: StateShape distr -> HMMCSV.CSVParser distr+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@@ -390,21 +459,18 @@    symbolFromCell = maybeRead  -instance-   (Shape.C sh, Class.Real prob, Show prob, Read prob, CSVSymbol symbol) =>-      ToCSV (Discrete symbol sh prob) where+instance (CSVSymbol symbol) => ToCSV (Discrete symbol) where    toCells (Discrete m) =       map          (\(symbol, probs) ->             cellFromSymbol symbol : HMMCSV.cellsFromVector probs) $-      Map.toAscList m+      Array.toAssociations $ Matrix.toRowArray m -instance-   (Shape.C sh, Class.Real prob, Show prob, Read prob, CSVSymbol symbol) =>-      FromCSV (Discrete symbol sh prob) where+instance (CSVSymbol symbol) => FromCSV (Discrete symbol) where    parseCells n =-      fmap (Discrete . Map.fromList) $-      HMMCSV.manyRowsUntilEnd $ parseSymbolProb 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) =>@@ -426,36 +492,30 @@                 return $ StorableArray.reshape sh v)  -instance-   (Shape.Indexed emiSh, Shape.Indexed stateSh,-    Class.Real a, Eq a, Show a, Read a) =>-      ToCSV (Gaussian emiSh stateSh a) where+instance (Shape.Indexed emiSh) => ToCSV (Gaussian emiSh) where    toCells (Gaussian params) =       List.intercalate [[]] $       map-         (\(center, covarianceChol, _) ->+         (\(_, center, covarianceChol) ->             HMMCSV.cellsFromVector center :             HMMCSV.cellsFromSquare (Triangular.toSquare covarianceChol)) $       Array.toList params -instance-   (emiSh ~ Matrix.ZeroInt, Shape.Indexed stateSh,-    Class.Real a, Eq a, Show a, Read a) =>-      FromCSV (Gaussian emiSh stateSh a) where+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 . fst3) gs+      let sizes = map (vectorDim . snd3) gs       HMMCSV.assert (ListHT.allEqual sizes) $          printf "dimensions of emissions mismatch: %s" (show sizes)-      return $ consGaussian $ Array.fromList sh gs+      return $ Gaussian $ Array.fromList sh gs  parseSingleGaussian ::-   (emiSh ~ Matrix.ZeroInt, Class.Real prob, Eq prob, Read prob) =>-   HMMCSV.CSVParser (Vector emiSh prob, UpperTriangular emiSh prob, prob)+   (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 <-@@ -475,6 +535,8 @@ -} isUpperTriang ::    (Shape.C sh, Class.Real a, Eq a) =>-   SquareMatrix sh a -> UpperTriangular sh a -> Bool+   Matrix.Square sh a -> Triangular.Upper sh a -> Bool isUpperTriang m mt =-   Vector.toList m == Vector.toList (Triangular.toSquare mt)+   Vector.toList (ArrMatrix.toVector m)+   ==+   Vector.toList (ArrMatrix.toVector (Triangular.toSquare mt))

@@ -46,7 +46,7 @@ hmm :: HMM hmm =    HMM.Cons {-      HMM.initial = normalizeProb $ Vector.constant stateSet 1,+      HMM.initial = normalizeProb $ Vector.one stateSet,       HMM.transition =          squareFromLists stateSet $             stateVector 0.9 0.0 0.0 0.1 :

@@ -40,7 +40,7 @@ hmm :: HMM hmm =    HMM.Cons {-      HMM.initial = normalizeProb $ Vector.constant stateSet 1,+      HMM.initial = normalizeProb $ Vector.one stateSet,       HMM.transition =          squareFromLists stateSet $             stateVector 0.9 0.0 0.0 0.1 :

@@ -7,17 +7,14 @@  import qualified Numeric.LAPACK.Vector as Vector import Numeric.LAPACK.Vector (Vector)-import Numeric.LAPACK.Format (Format(format))  import qualified Data.Array.Comfort.Shape as Shape -import qualified Text.PrettyPrint.Boxes as TextBox import Text.Read.HT (maybeRead)  import Control.DeepSeq (NFData(rnf)) import Control.Monad (liftM2) -import qualified Data.Map as Map import qualified Data.NonEmpty as NonEmpty import qualified Data.List.HT as ListHT import Data.NonEmpty ((!:))@@ -31,9 +28,6 @@    rnf Red = ()    rnf _ = () -instance Format Color where-   format _fmt = TextBox.text . show- {- | Using 'show' and 'read' is not always a good choice since they must format and parse Haskell expressions@@ -67,9 +61,9 @@             stateVector 0.0 0.0 0.2 0.8 :             [],       HMM.distribution =-         Distr.Discrete $ Map.fromList $-            (Red,    stateVector 1 0 0 0) :-            (Yellow, stateVector 0 1 0 1):+         Distr.discreteFromList $+            (Red,    stateVector 1 0 0 0) !:+            (Yellow, stateVector 0 1 0 1) :             (Green,  stateVector 0 0 1 0) :             []    }@@ -86,8 +80,8 @@             stateVector 0.2 0.2 0.3 0.3 :             [],       HMM.distribution =-         Distr.Discrete $ Map.fromList $-            (Red,    stateVector 0.6 0.2 0.2 0.2) :+         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) :             []

@@ -1,6 +1,4 @@ {-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE UndecidableInstances #-} module Math.HiddenMarkovModel.Named (    T(..),    Discrete,@@ -42,30 +40,31 @@ Although 'nameFromStateMap' and 'stateFromNameMap' are exported you must be careful to keep them consistent when you alter them. -}-data T distr sh ix prob =+data T typ sh ix prob =    Cons {-      model :: HMM.T distr sh prob,+      model :: HMM.T typ sh prob,       nameFromStateMap :: Array sh String,       stateFromNameMap :: Map String ix    }    deriving (Show) +type Simple typ sh prob = T typ sh (Shape.Index sh) prob type Discrete symbol stateSh prob =-      T (Distr.Discrete symbol stateSh prob) stateSh (Shape.Index stateSh) prob+      Simple (Distr.Discrete symbol) stateSh prob type Gaussian emiSh stateSh a =-      T (Distr.Gaussian emiSh stateSh a) stateSh (Shape.Index stateSh) a+      Simple (Distr.Gaussian emiSh) stateSh a   instance-   (NFData distr, NFData sh, NFData ix, NFData prob,+   (Distr.NFData typ, NFData sh, NFData ix, NFData prob,     Shape.C sh, Storable prob) =>-      NFData (T distr sh ix prob) where+      NFData (T typ sh ix prob) where    rnf hmm = rnf (model hmm, nameFromStateMap hmm, stateFromNameMap hmm)   fromModelAndNames ::-   (Shape.Indexed sh, Shape.Index sh ~ state) =>-   HMM.T distr sh prob -> [String] -> T distr sh state prob+   (Shape.Indexed sh) =>+   HMM.T typ sh prob -> [String] -> Simple typ sh prob fromModelAndNames md names =    let m = Array.fromList (StorableArray.shape $ HMM.initial md) names    in  Cons {@@ -82,26 +81,24 @@   toCSV ::-   (Distr.ToCSV distr, Shape.Indexed sh, Class.Real prob, Show prob) =>-   T distr sh ix prob -> String+   (Distr.ToCSV typ, Shape.Indexed sh, Class.Real prob, Show prob) =>+   Simple typ sh prob -> String toCSV hmm =    CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $       Array.toList (nameFromStateMap hmm) : HMM.toCells (model hmm)  fromCSV ::-   (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh,-    Shape.Indexed stateSh, Shape.Index stateSh ~ state,+   (Distr.FromCSV typ, Shape.Indexed stateSh, Eq stateSh,     Class.Real prob, Read prob) =>    (Int -> stateSh) ->-   String -> ME.Exceptional String (T distr stateSh state prob)+   String -> ME.Exceptional String (Simple typ stateSh prob) fromCSV makeShape =    MS.evalStateT (parseCSV makeShape) . map HMMCSV.fixShortRow . CSV.parseCSV  parseCSV ::-   (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh,-    Shape.Indexed stateSh, Shape.Index stateSh ~ state,+   (Distr.FromCSV typ, Shape.Indexed stateSh, Eq stateSh,     Class.Real prob, Read prob) =>-   (Int -> stateSh) -> HMMCSV.CSVParser (T distr stateSh state prob)+   (Int -> stateSh) -> HMMCSV.CSVParser (Simple typ stateSh prob) parseCSV makeShape = do    names <- HMMCSV.parseStringList =<< HMMCSV.getRow    let duplicateNames =

@@ -10,12 +10,12 @@ import qualified Math.HiddenMarkovModel.Distribution as Distr import Math.HiddenMarkovModel.Private           (T(..), Trained(..), emission,-           biscaleTransition, matrixMaxMul, sumTransitions)-import Math.HiddenMarkovModel.Utility-         (SquareMatrix, normalizeFactor, normalizeProb)+           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.Matrix ((-*#), (#*|)) import Numeric.LAPACK.Vector (Vector)  import qualified Numeric.Netlib.Class as Class@@ -23,15 +23,12 @@ import qualified Control.Functor.HT as Functor  import qualified Data.Array.Comfort.Storable as StorableArray-import qualified Data.Array.Comfort.Boxed as Array import qualified Data.Array.Comfort.Shape as Shape  import qualified Data.NonEmpty.Class as NonEmptyC import qualified Data.NonEmpty as NonEmpty import qualified Data.Foldable as Fold-import qualified Data.List as List-import Data.Traversable (Traversable, mapAccumL)-import Data.Tuple.HT (mapFst, mapSnd, swap)+import Data.Traversable (Traversable)   {- |@@ -40,43 +37,43 @@ that it may be rounded to zero in the choosen number type. -} logLikelihood ::-   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh, Floating prob,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Floating prob,+    Class.Real prob, Distr.Emission typ prob ~ emission,     Traversable f) =>-   T distr sh prob -> NonEmpty.T f emission -> prob+   T typ sh prob -> NonEmpty.T f emission -> prob logLikelihood hmm = Fold.sum . fmap (log . fst) . alpha hmm  alpha ::-   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+   (Distr.EmissionProb typ, Shape.C sh, Eq sh,+    Class.Real prob, Distr.Emission typ prob ~ emission,     Traversable f) =>-   T distr sh prob ->+   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))+          (\(_,alphai) xi -> normMulEmiss xi (transition hmm #*| alphai))           (normMulEmiss x (initial hmm))           xs  beta ::-   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+   (Distr.EmissionProb typ, Shape.C sh, Eq sh,+    Class.Real prob, Distr.Emission typ prob ~ emission,     Traversable f, NonEmptyC.Reverse f) =>-   T distr sh prob ->+   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.constant (StorableArray.shape $ initial hmm) 1)+         Vector.mul (emission hmm xi) betai -*# transition hmm)+      (Vector.one $ StorableArray.shape $ initial hmm)  alphaBeta ::-   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+   (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 distr sh prob ->+   T typ sh prob ->    NonEmpty.T f emission ->    (NonEmpty.T f (prob, Vector sh prob), NonEmpty.T f (Vector sh prob)) alphaBeta hmm xs =@@ -86,19 +83,19 @@   xiFromAlphaBeta ::-   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,+   (Distr.EmissionProb typ, Shape.C sh, Eq sh,+    Class.Real prob, Distr.Emission typ prob ~ emission,     Traversable f, NonEmptyC.Zip f) =>-   T distr sh prob ->+   T typ sh prob ->    NonEmpty.T f emission ->    NonEmpty.T f (prob, Vector sh prob) ->    NonEmpty.T f (Vector sh prob) ->-   f (SquareMatrix 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 ->-             Vector.scale (recip c1) $ biscaleTransition hmm x alpha0 beta1)+             Matrix.scale (recip c1) $ biscaleTransition hmm x alpha0 beta1)           (NonEmpty.tail xs)           (NonEmpty.init alphas)           (NonEmpty.tail cs)@@ -119,22 +116,10 @@ It is found using the Viterbi algorithm. -} reveal ::-   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh,-    Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ state,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,-    Traversable f, NonEmptyC.Reverse f) =>-   T distr sh prob -> NonEmpty.T f emission -> NonEmpty.T f state-reveal hmm (NonEmpty.Cons x xs) =-   fmap (Shape.revealIndex (StorableArray.shape $ initial hmm)) $-   uncurry (NonEmpty.scanr (StorableArray.!)) $-   mapFst-      (fst . Vector.argAbsMaximum .-       StorableArray.mapShape Shape.Deferred) $-   mapAccumL-      (\alphai xi ->-         swap $ mapSnd (Vector.mul (emission hmm xi)) $-         matrixMaxMul (transition hmm) $ normalizeProb alphai)-      (Vector.mul (emission hmm x) (initial hmm)) xs+   (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   {- |@@ -155,9 +140,9 @@ This is done by the Baum-Welch algorithm. -} trainUnsupervised ::-   (Distr.Estimate tdistr distr, Distr.StateShape distr ~ sh, Eq sh,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>-   T distr sh prob -> NonEmpty.T [] emission -> Trained tdistr sh prob+   (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@@ -168,8 +153,5 @@           trainedTransition =              sumTransitions hmm $ xiFromAlphaBeta hmm xs alphas betas,           trainedDistribution =-             Distr.accumulateEmissions $-             Array.fromList (StorableArray.shape zeta0) $-             map (zip (NonEmpty.flatten xs)) $-             List.transpose $ map Vector.toList $ NonEmpty.flatten zetas+             Distr.accumulateEmissionVectors $ NonEmptyC.zip xs zetas        }

@@ -30,8 +30,10 @@ import qualified Math.HiddenMarkovModel.Distribution as Distr import qualified Math.HiddenMarkovModel as HMM import Math.HiddenMarkovModel.Private (Trained(..))-import Math.HiddenMarkovModel.Utility (SquareMatrix, squareConstant)+import Math.HiddenMarkovModel.Utility (squareConstant) +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 @@ -44,16 +46,16 @@ import Data.FixedLength ((!:))  import qualified Type.Data.Num.Unary.Literal as TypeNum-import Type.Base.Proxy (Proxy(Proxy)) -import qualified Data.Map as Map+import qualified Data.NonEmpty.Map as NonEmptyMap+import qualified Data.NonEmpty as NonEmpty import Data.Semigroup (Semigroup, (<>), stimes)  import Prelude hiding (replicate)   newtype T sh prob =-   Cons (sh -> (Shape.Index sh, SquareMatrix sh prob, Shape.Index sh))+   Cons (sh -> (Shape.Index sh, Matrix.Square sh prob, Shape.Index sh))  atom ::    (Shape.Indexed sh, Shape.Index sh ~ state, Class.Real prob) =>@@ -77,7 +79,7 @@    Cons $ \n ->       case (f n, g n) of          ((sai, ma, sao), (sbi, mb, sbo)) ->-            (sai, increment (sbi,sao) 1 $ Vector.add ma mb, sbo)+            (sai, increment (sbi,sao) 1 $ Matrix.add ma mb, sbo)  replicate ::    (Shape.Indexed sh, Class.Real prob) => Int -> T sh prob -> T sh prob@@ -86,25 +88,26 @@       case f sh of          (si, m, so) ->             let k = fromIntegral ki-            in  (si, increment (si,so) (k-1) $ Vector.scale k m, so)+            in  (si, increment (si,so) (k-1) $ Matrix.scale k m, so)  increment ::    (Shape.Indexed sh, Shape.Index sh ~ state, Class.Real a) =>-   (state, state) -> a -> SquareMatrix sh a -> SquareMatrix sh a-increment (i,j) x m  =  StorableArray.accumulate (+) m [((i,j), x)]+   (state, state) -> a -> Matrix.Square sh a -> Matrix.Square sh a+increment (i,j) x =+   ArrMatrix.lift1 $ flip (StorableArray.accumulate (+)) [((i,j), x)]   finish ::-   (Shape.Indexed sh, Class.Real prob) =>-   sh -> tdistr -> T sh prob -> Trained tdistr sh prob-finish sh tdistr (Cons f) =-   case f sh of-      (si, m, _so) ->-         Trained {-            trainedInitial = StorableArray.fromAssociations sh 0 [(si,1)],-            trainedTransition = m,-            trainedDistribution = tdistr-         }+   (Distr.Info typ, Shape.Indexed sh, Class.Real prob) =>+   Distr.Trained typ sh prob -> T sh prob -> Trained typ sh prob+finish tdistr (Cons f) =+   let sh = Distr.statesShapeTrained tdistr+       (si, m, _so) = f sh+   in Trained {+         trainedInitial = Vector.unit sh si,+         trainedTransition = m,+         trainedDistribution = tdistr+      }   _example :: HMM.DiscreteTrained Char (ShapeStatic.ZeroBased TypeNum.U2) Double@@ -112,10 +115,9 @@    let a = atom FL.i0        b = atom FL.i1        distr =-          Distr.DiscreteTrained $ Map.fromList $-          ('a', ShapeStatic.vector $ 1!:2!:FL.end) :+          Distr.DiscreteTrained $ NonEmptyMap.fromList $+          ('a', ShapeStatic.vector $ 1!:2!:FL.end) NonEmpty.!:           ('b', ShapeStatic.vector $ 4!:3!:FL.end) :           ('c', ShapeStatic.vector $ 0!:1!:FL.end) :           []-   in finish (ShapeStatic.ZeroBased Proxy) distr $-      replicate 5 $ replicate 10 a <> replicate 20 b+   in finish distr $ replicate 5 $ replicate 10 a <> replicate 20 b

@@ -1,15 +1,16 @@ {-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE UndecidableInstances #-} module Math.HiddenMarkovModel.Private where  import qualified Math.HiddenMarkovModel.Distribution as Distr import qualified Math.HiddenMarkovModel.CSV as HMMCSV-import Math.HiddenMarkovModel.Utility (SquareMatrix, diagonal)+import 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.Matrix ((-*#), (##*#), (#*##), (#*|)) import Numeric.LAPACK.Vector (Vector)  import qualified Numeric.Netlib.Class as Class@@ -20,15 +21,15 @@ import Foreign.Storable (Storable)  import qualified Data.Array.Comfort.Storable as StorableArray-import qualified Data.Array.Comfort.Boxed as Array import qualified Data.Array.Comfort.Shape as Shape  import qualified Data.NonEmpty.Class as NonEmptyC import qualified Data.NonEmpty as NonEmpty import qualified Data.Semigroup as Sg import qualified Data.List as List+import Data.Semigroup ((<>)) import Data.Traversable (Traversable, mapAccumL)-import Data.Tuple.HT (mapPair, mapFst, mapSnd, swap)+import Data.Tuple.HT (mapFst, mapSnd, swap)   {- |@@ -46,88 +47,84 @@ with respect to popular HMM descriptions. But I think this is the natural orientation, because this way you can write \"transition matrix times probability column vector\".--The type has two type parameters,-although the one for the distribution would be enough.-However, replacing @prob@ by @Distr.Probability distr@-would prohibit the derived Show and Read instances. -}-data T distr sh prob =+data T typ sh prob =    Cons {       initial :: Vector sh prob,-      transition :: SquareMatrix sh prob,-      distribution :: distr+      transition :: Matrix.Square sh prob,+      distribution :: Distr.T typ sh prob    }    deriving (Show)  instance-   (NFData distr, NFData sh, NFData prob, Storable prob) =>-      NFData (T distr sh prob) where+   (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-   (Class.Real prob, Format.FormatArray sh, Format.Format distr) =>-      Format.Format (T distr sh prob) where+   (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 ::-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>-   T distr sh prob -> emission -> Vector sh prob+   (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 ::-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,-    Traversable f) =>-   T distr sh prob -> NonEmpty.T f emission -> prob+   (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 ::-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,-    Traversable f) =>-   T distr sh prob ->-   NonEmpty.T f emission -> NonEmpty.T f (Vector sh prob)+   (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))+      (\alphai xi -> Vector.mul (emission hmm xi) (transition hmm #*| alphai))       (Vector.mul (emission hmm x) (initial hmm))       xs   backward ::-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,-    Traversable f) =>-   T distr sh prob -> NonEmpty.T f emission -> prob+   (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.sum $-   Vector.mul (initial hmm) $+   Vector.dot (initial hmm) $    Vector.mul (emission hmm x) $    NonEmpty.head $ beta hmm xs  beta ::-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,-    Traversable f) =>-   T distr sh prob ->-   f emission -> NonEmpty.T f (Vector sh prob)+   (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.constant (StorableArray.shape $ initial hmm) 1)+      (\xi betai -> Vector.mul (emission hmm xi) betai -*# transition hmm)+      (Vector.one $ StorableArray.shape $ initial hmm)   alphaBeta ::-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,-    Traversable f) =>-   T distr sh prob ->+   (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 =@@ -139,29 +136,28 @@   biscaleTransition ::-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr,-    Distr.EmissionProb distr, Distr.Probability distr ~ prob) =>-   T distr sh prob -> Distr.Emission distr ->-   Vector sh prob -> Vector sh prob -> SquareMatrix sh prob+   (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 (Vector.mul (emission hmm x) beta1)+    #*##+    transition hmm)+   ##*#    diagonal alpha0  xiFromAlphaBeta ::-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>-   T distr sh prob -> prob ->+   (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) ->-   [SquareMatrix sh prob]+   [Matrix.Square sh prob] xiFromAlphaBeta hmm recipLikelihood xs alphas betas =    zipWith3       (\x alpha0 beta1 ->-         Vector.scale recipLikelihood $+         Matrix.scale recipLikelihood $          biscaleTransition hmm x alpha0 beta1)       (NonEmpty.tail xs)       (NonEmpty.init alphas)@@ -169,7 +165,7 @@  zetaFromXi ::    (Shape.C sh, Eq sh, Class.Real prob) =>-   [SquareMatrix sh prob] -> [Vector sh prob]+   [Matrix.Square sh prob] -> [Vector sh prob] zetaFromXi = map Matrix.columnSums  zetaFromAlphaBeta ::@@ -191,35 +187,41 @@ than the smallest representable number. -} reveal ::-   (Shape.InvIndexed sh, Shape.Index sh ~ state,-    Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,-    Traversable f) =>-   T distr sh prob -> NonEmpty.T f emission -> NonEmpty.T f state-reveal hmm (NonEmpty.Cons x xs) =-   fmap (Shape.revealIndex (StorableArray.shape $ initial hmm)) $+   (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 .-       StorableArray.mapShape Shape.Deferred) $+   mapFst (fst . Vector.argAbsMaximum) $    mapAccumL       (\alphai xi ->          swap $ mapSnd (Vector.mul (emission hmm xi)) $-         matrixMaxMul (transition hmm) alphai)+         matrixMaxMul (transition hmm) $ normalize alphai)       (Vector.mul (emission hmm x) (initial hmm)) xs  matrixMaxMul ::-   (Shape.Indexed sh, Eq sh, Shape.Index sh ~ ix, Class.Real a) =>-   SquareMatrix sh a -> Vector sh a ->-   (Vector (Shape.Deferred sh) (Shape.DeferredIndex ix), Vector sh a)-matrixMaxMul m v =-   let sh = StorableArray.shape v-   in mapPair (Vector.fromList (Shape.Deferred sh), Vector.fromList sh) $-      unzip $-      map (Vector.argAbsMaximum .-           StorableArray.mapShape Shape.Deferred .-           Vector.mul v) $-      Matrix.toRows m+   (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   @@ -239,35 +241,35 @@  * derive it from state sequence patterns, cf. "Math.HiddenMarkovModel.Pattern". -}-data Trained distr sh prob =+data Trained typ sh prob =    Trained {       trainedInitial :: Vector sh prob,-      trainedTransition :: SquareMatrix sh prob,-      trainedDistribution :: distr+      trainedTransition :: Matrix.Square sh prob,+      trainedDistribution :: Distr.Trained typ sh prob    }    deriving (Show)  instance-   (NFData distr, NFData sh, NFData prob, Storable prob) =>-      NFData (Trained distr sh prob) where+   (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 distr sh e -> [SquareMatrix sh e] -> SquareMatrix sh e+   T typ sh e -> [Matrix.Square sh e] -> Matrix.Square sh e sumTransitions hmm =-   List.foldl' Vector.add-      (Vector.constant (StorableArray.shape $ transition hmm) 0)+   List.foldl' Matrix.add $+   Matrix.zero $ ArrMatrix.shape $ transition hmm  {- | Baum-Welch algorithm -} trainUnsupervised ::-   (Distr.Estimate tdistr distr, Distr.StateShape distr ~ sh, Eq sh,-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>-   T distr sh prob -> NonEmpty.T [] emission -> Trained tdistr sh prob+   (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@@ -279,37 +281,31 @@              sumTransitions hmm $              xiFromAlphaBeta hmm recipLikelihood xs alphas betas,           trainedDistribution =-             Distr.accumulateEmissions $-             Array.fromList (StorableArray.shape zeta0) $-             map (zip (NonEmpty.flatten xs)) $-             List.transpose $ map Vector.toList $ NonEmpty.flatten zetas+             Distr.accumulateEmissionVectors $ NonEmptyC.zip xs zetas        }   mergeTrained ::-   (Shape.C sh, Eq sh,-    Distr.Estimate tdistr distr, Distr.Probability distr ~ prob) =>-   Trained tdistr sh prob -> Trained tdistr sh prob -> Trained tdistr sh prob+   (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 =-         Vector.add (trainedTransition hmm0) (trainedTransition hmm1),+         Matrix.add (trainedTransition hmm0) (trainedTransition hmm1),       trainedDistribution =-         Distr.combine-            (trainedDistribution hmm0) (trainedDistribution hmm1)+         trainedDistribution hmm0 <> trainedDistribution hmm1    }  instance-   (Shape.C sh, Eq sh,-    Distr.Estimate tdistr distr, Distr.Probability distr ~ prob) =>-      Sg.Semigroup (Trained tdistr sh prob) where+   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>+      Sg.Semigroup (Trained typ sh prob) where    (<>) = mergeTrained   toCells ::-   (Distr.ToCSV distr, Shape.Indexed sh, Class.Real prob, Show prob) =>-   T distr sh prob -> [[String]]+   (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) ++@@ -317,9 +313,9 @@    (Distr.toCells $ distribution hmm)  parseCSV ::-   (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh, Shape.C stateSh,+   (Distr.FromCSV typ, Shape.C stateSh, Eq stateSh,     Class.Real prob, Read prob) =>-   (Int -> stateSh) -> HMMCSV.CSVParser (T distr stateSh prob)+   (Int -> stateSh) -> HMMCSV.CSVParser (T typ stateSh prob) parseCSV makeShape = do    v <-       StorableArray.mapShape (makeShape . Shape.zeroBasedSize) <$>

@@ -11,8 +11,10 @@ import qualified Math.HiddenMarkovModel.Normalized as Normalized import qualified Math.HiddenMarkovModel.Private as Priv import qualified Math.HiddenMarkovModel.Distribution as Distr-import Math.HiddenMarkovModel.Utility (SquareMatrix, squareFromLists, distance)+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)@@ -30,10 +32,12 @@ 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)@@ -53,8 +57,8 @@             stateVector 0.1 0.1 0.2 0.8 :             [],       HMM.distribution =-         Distr.Discrete $ Map.fromList $-            ('a', stateVector 1 0 0 0) :+         Distr.discreteFromList $+            ('a', stateVector 1 0 0 0) !:             ('b', stateVector 0 1 0 1) :             ('c', stateVector 0 0 1 0) :             []@@ -83,8 +87,7 @@             1 -> True             _ -> error "invalid emission probability (must be 0 or 1)") $    Vector.toList $-   Map.findWithDefault (error "invalid character") c $-   case HMM.distribution hmm of Distr.Discrete m -> m+   case HMM.distribution hmm of Distr.Discrete m -> Matrix.takeRow m c  {- | Should all be equal.@@ -140,7 +143,7 @@ {- | Lists should be equal -}-xis :: ([SquareMatrix StateSet Double], [SquareMatrix StateSet Double])+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,@@ -153,7 +156,8 @@ xisDiff :: (Bool, Double) xisDiff =    case xis of-      (x0,x1) -> (length x0 == length x1, maximum $ zipWith distance x0 x1)+      (x0,x1) ->+         (length x0 == length x1, maximum $ zipWith matrixDistance x0 x1)   reveal :: Bool@@ -172,13 +176,16 @@ trainUnsupervisedDiff =    case trainUnsupervised of       (hmm0,hmm1) ->-         (distance (Priv.trainedTransition hmm0) (Priv.trainedTransition 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) ->-                (Map.size m0 == Map.size m1,-                 Fold.maximum $ Map.intersectionWith distance m0 m1))+                (NonEmptyMap.size m0 == NonEmptyMap.size m1,+                 Fold.maximum $+                 Map.intersectionWith distance+                    (NonEmptyMap.flatten m0) (NonEmptyMap.flatten m1)))   nonEmptyScanr :: Int -> [Int] -> Bool

@@ -4,10 +4,12 @@ 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.Vector (Vector)+import Numeric.LAPACK.Matrix.Array (ArrayMatrix)+import Numeric.LAPACK.Vector (Vector, (.*|))  import qualified Numeric.Netlib.Class as Class @@ -22,15 +24,13 @@ import qualified Control.Monad.Trans.State as MS  -type SquareMatrix sh = Square.Square sh- normalizeProb :: (Shape.C sh, Class.Real a) => Vector sh a -> Vector sh a normalizeProb = snd . normalizeFactor  normalizeFactor :: (Shape.C sh, Class.Real a) => Vector sh a -> (a, Vector sh a) normalizeFactor xs =    let c = Vector.sum xs-   in  (c, Vector.scale (recip c) xs)+   in  (c, recip c .*| xs)  -- see htam:Stochastic randomItemProp ::@@ -57,11 +57,13 @@   squareConstant ::-   (Shape.C sh, Class.Real a) => sh -> a -> SquareMatrix sh a-squareConstant = Vector.constant . MatrixShape.square MatrixShape.RowMajor+   (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] -> SquareMatrix sh a+   (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 @@ -79,3 +81,8 @@    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)