module ELynx.Data.MarkovProcess.RateMatrix
( RateMatrix
, ExchangeabilityMatrix
, StationaryDistribution
, totalRate
, normalize
, normalizeWith
, setDiagonal
, toExchangeabilityMatrix
, fromExchangeabilityMatrix
, getStationaryDistribution
) where
import Numeric.LinearAlgebra hiding (normalize)
import Prelude hiding ((<>))
import ELynx.Tools.Equality
import ELynx.Tools.LinearAlgebra
import ELynx.Tools.Vector
type RateMatrix = Matrix R
type ExchangeabilityMatrix = Matrix R
type StationaryDistribution = Vector R
totalRate :: StationaryDistribution -> RateMatrix -> Double
totalRate d m = norm_1 $ d <# matrixSetDiagToZero m
normalize :: RateMatrix -> RateMatrix
normalize m = normalizeWith (getStationaryDistribution m) m
normalizeWith :: StationaryDistribution -> RateMatrix -> RateMatrix
normalizeWith d m = scale (1.0 / totalRate d m) m
setDiagonal :: RateMatrix -> RateMatrix
setDiagonal m = diagZeroes - diag (fromList rowSums)
where diagZeroes = matrixSetDiagToZero m
rowSums = map norm_1 $ toRows diagZeroes
toExchangeabilityMatrix :: RateMatrix -> StationaryDistribution -> ExchangeabilityMatrix
toExchangeabilityMatrix m f = m <> diag oneOverF
where oneOverF = cmap (1.0/) f
fromExchangeabilityMatrix :: ExchangeabilityMatrix -> StationaryDistribution -> RateMatrix
fromExchangeabilityMatrix em d = setDiagonal $ em <> diag d
getStationaryDistribution :: RateMatrix -> StationaryDistribution
getStationaryDistribution m =
if magnitude (eVals ! i) `nearlyEq` 0
then normalizeSumVec 1.0 distReal
else error "Could not retrieve stationary distribution."
where
(eVals, eVecs) = eig (tr m)
i = minIndex eVals
distComplex = toColumns eVecs !! i
distReal = cmap realPart distComplex