{-# LANGUAGE FlexibleContexts #-}
module ELynx.MarkovProcess.RateMatrix
( RateMatrix,
ExchangeabilityMatrix,
StationaryDistribution,
isValid,
normalizeSD,
totalRate,
totalRateWith,
normalize,
normalizeWith,
setDiagonal,
toExchangeabilityMatrix,
fromExchangeabilityMatrix,
getStationaryDistribution,
exchFromListLower,
exchFromListUpper,
)
where
import qualified Data.Vector.Storable as V
import Numeric.LinearAlgebra hiding (normalize)
import Numeric.SpecFunctions
import Prelude hiding ((<>))
type RateMatrix = Matrix R
type ExchangeabilityMatrix = Matrix R
type StationaryDistribution = Vector R
epsRelaxed :: Double
epsRelaxed :: Double
epsRelaxed = Double
1e-5
isValid :: StationaryDistribution -> Bool
isValid :: Vector Double -> Bool
isValid Vector Double
d = Double
epsRelaxed forall a. Ord a => a -> a -> Bool
> forall a. Num a => a -> a
abs (forall a. Normed a => a -> Double
norm_1 Vector Double
d forall a. Num a => a -> a -> a
- Double
1.0)
normalizeSD :: StationaryDistribution -> StationaryDistribution
normalizeSD :: Vector Double -> Vector Double
normalizeSD Vector Double
d = Vector Double
d forall a. Fractional a => a -> a -> a
/ forall (c :: * -> *) e. Container c e => e -> c e
scalar (forall a. Normed a => a -> Double
norm_1 Vector Double
d)
matrixSetDiagToZero :: Matrix R -> Matrix R
matrixSetDiagToZero :: Matrix Double -> Matrix Double
matrixSetDiagToZero Matrix Double
m = Matrix Double
m forall a. Num a => a -> a -> a
- forall a. (Num a, Element a) => Vector a -> Matrix a
diag (forall t. Element t => Matrix t -> Vector t
takeDiag Matrix Double
m)
{-# INLINE matrixSetDiagToZero #-}
totalRateWith :: StationaryDistribution -> RateMatrix -> Double
totalRateWith :: Vector Double -> Matrix Double -> Double
totalRateWith Vector Double
d Matrix Double
m = forall a. Normed a => a -> Double
norm_1 forall a b. (a -> b) -> a -> b
$ Vector Double
d forall t. Numeric t => Vector t -> Matrix t -> Vector t
<# Matrix Double -> Matrix Double
matrixSetDiagToZero Matrix Double
m
totalRate :: RateMatrix -> Double
totalRate :: Matrix Double -> Double
totalRate Matrix Double
m = Vector Double -> Matrix Double -> Double
totalRateWith (Matrix Double -> Vector Double
getStationaryDistribution Matrix Double
m) Matrix Double
m
normalize :: RateMatrix -> RateMatrix
normalize :: Matrix Double -> Matrix Double
normalize Matrix Double
m = Vector Double -> Matrix Double -> Matrix Double
normalizeWith (Matrix Double -> Vector Double
getStationaryDistribution Matrix Double
m) Matrix Double
m
normalizeWith :: StationaryDistribution -> RateMatrix -> RateMatrix
normalizeWith :: Vector Double -> Matrix Double -> Matrix Double
normalizeWith Vector Double
d Matrix Double
m = forall t (c :: * -> *). Linear t c => t -> c t -> c t
scale (Double
1.0 forall a. Fractional a => a -> a -> a
/ Vector Double -> Matrix Double -> Double
totalRateWith Vector Double
d Matrix Double
m) Matrix Double
m
setDiagonal :: RateMatrix -> RateMatrix
setDiagonal :: Matrix Double -> Matrix Double
setDiagonal Matrix Double
m = Matrix Double
diagZeroes forall a. Num a => a -> a -> a
- forall a. (Num a, Element a) => Vector a -> Matrix a
diag (forall a. Storable a => [a] -> Vector a
fromList [Double]
rowSums)
where
diagZeroes :: Matrix Double
diagZeroes = Matrix Double -> Matrix Double
matrixSetDiagToZero Matrix Double
m
rowSums :: [Double]
rowSums = forall a b. (a -> b) -> [a] -> [b]
map forall a. Normed a => a -> Double
norm_1 forall a b. (a -> b) -> a -> b
$ forall t. Element t => Matrix t -> [Vector t]
toRows Matrix Double
diagZeroes
toExchangeabilityMatrix ::
RateMatrix -> StationaryDistribution -> ExchangeabilityMatrix
toExchangeabilityMatrix :: Matrix Double -> Vector Double -> Matrix Double
toExchangeabilityMatrix Matrix Double
m Vector Double
f = Matrix Double
m forall t. Numeric t => Matrix t -> Matrix t -> Matrix t
<> forall a. (Num a, Element a) => Vector a -> Matrix a
diag Vector Double
oneOverF
where
oneOverF :: Vector Double
oneOverF = forall b (c :: * -> *) e.
(Element b, Container c e) =>
(e -> b) -> c e -> c b
cmap (Double
1.0 forall a. Fractional a => a -> a -> a
/) Vector Double
f
fromExchangeabilityMatrix ::
ExchangeabilityMatrix -> StationaryDistribution -> RateMatrix
fromExchangeabilityMatrix :: Matrix Double -> Vector Double -> Matrix Double
fromExchangeabilityMatrix Matrix Double
em Vector Double
d = Matrix Double -> Matrix Double
setDiagonal forall a b. (a -> b) -> a -> b
$ Matrix Double
em forall t. Numeric t => Matrix t -> Matrix t -> Matrix t
<> forall a. (Num a, Element a) => Vector a -> Matrix a
diag Vector Double
d
eps :: Double
eps :: Double
eps = Double
1e-12
normalizeSumVec :: V.Vector Double -> V.Vector Double
normalizeSumVec :: Vector Double -> Vector Double
normalizeSumVec Vector Double
v = forall a b.
(Storable a, Storable b) =>
(a -> b) -> Vector a -> Vector b
V.map (forall a. Fractional a => a -> a -> a
/ Double
s) Vector Double
v
where
s :: Double
s = forall a. (Storable a, Num a) => Vector a -> a
V.sum Vector Double
v
{-# INLINE normalizeSumVec #-}
getStationaryDistribution :: RateMatrix -> StationaryDistribution
getStationaryDistribution :: Matrix Double -> Vector Double
getStationaryDistribution Matrix Double
m =
if Double
eps forall a. Ord a => a -> a -> Bool
> forall a. Num a => a -> a
abs (forall a. RealFloat a => Complex a -> a
magnitude (Vector (Complex Double)
eVals forall c t. Indexable c t => c -> Int -> t
! Int
i))
then Vector Double -> Vector Double
normalizeSumVec Vector Double
distReal
else forall a. HasCallStack => [Char] -> a
error [Char]
"getStationaryDistribution: Could not retrieve stationary distribution."
where
(Vector (Complex Double)
eVals, Matrix (Complex Double)
eVecs) = forall t.
Field t =>
Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
eig (forall m mt. Transposable m mt => m -> mt
tr Matrix Double
m)
i :: IndexOf Vector
i = forall (c :: * -> *) e. Container c e => c e -> IndexOf c
minIndex Vector (Complex Double)
eVals
distComplex :: Vector (Complex Double)
distComplex = forall t. Element t => Matrix t -> [Vector t]
toColumns Matrix (Complex Double)
eVecs forall a. [a] -> Int -> a
!! Int
i
distReal :: Vector Double
distReal = forall b (c :: * -> *) e.
(Element b, Container c e) =>
(e -> b) -> c e -> c b
cmap forall a. Complex a -> a
realPart Vector (Complex Double)
distComplex
ijToKLower :: Int -> Int -> Int
ijToKLower :: Int -> Int -> Int
ijToKLower Int
i Int
j
| Int
i forall a. Ord a => a -> a -> Bool
> Int
j = forall a b. (RealFrac a, Integral b) => a -> b
round (Int
i Int -> Int -> Double
`choose` Int
2) forall a. Num a => a -> a -> a
+ Int
j
| Bool
otherwise = forall a. HasCallStack => [Char] -> a
error [Char]
"ijToKLower: not defined for upper triangular matrix."
ijToKUpper :: Int -> Int -> Int -> Int
ijToKUpper :: Int -> Int -> Int -> Int
ijToKUpper Int
n Int
i Int
j
| Int
i forall a. Ord a => a -> a -> Bool
< Int
j = Int
i forall a. Num a => a -> a -> a
* (Int
n forall a. Num a => a -> a -> a
- Int
2) forall a. Num a => a -> a -> a
- forall a b. (RealFrac a, Integral b) => a -> b
round (Int
i Int -> Int -> Double
`choose` Int
2) forall a. Num a => a -> a -> a
+ Int
j forall a. Num a => a -> a -> a
- Int
1
| Bool
otherwise = forall a. HasCallStack => [Char] -> a
error [Char]
"ijToKUpper: not defined for lower triangular matrix."
fromListBuilderLower :: RealFrac a => [a] -> a -> a -> a
fromListBuilderLower :: forall a. RealFrac a => [a] -> a -> a -> a
fromListBuilderLower [a]
es a
i a
j
| a
i forall a. Ord a => a -> a -> Bool
> a
j = [a]
es forall a. [a] -> Int -> a
!! Int -> Int -> Int
ijToKLower Int
iI Int
jI
| a
i forall a. Eq a => a -> a -> Bool
== a
j = a
0.0
| a
i forall a. Ord a => a -> a -> Bool
< a
j = [a]
es forall a. [a] -> Int -> a
!! Int -> Int -> Int
ijToKLower Int
jI Int
iI
| Bool
otherwise =
forall a. HasCallStack => [Char] -> a
error
[Char]
"Float indices could not be compared during matrix creation."
where
iI :: Int
iI = forall a b. (RealFrac a, Integral b) => a -> b
round a
i :: Int
jI :: Int
jI = forall a b. (RealFrac a, Integral b) => a -> b
round a
j :: Int
fromListBuilderUpper :: RealFrac a => Int -> [a] -> a -> a -> a
fromListBuilderUpper :: forall a. RealFrac a => Int -> [a] -> a -> a -> a
fromListBuilderUpper Int
n [a]
es a
i a
j
| a
i forall a. Ord a => a -> a -> Bool
< a
j = [a]
es forall a. [a] -> Int -> a
!! Int -> Int -> Int -> Int
ijToKUpper Int
n Int
iI Int
jI
| a
i forall a. Eq a => a -> a -> Bool
== a
j = a
0.0
| a
i forall a. Ord a => a -> a -> Bool
> a
j = [a]
es forall a. [a] -> Int -> a
!! Int -> Int -> Int -> Int
ijToKUpper Int
n Int
jI Int
iI
| Bool
otherwise =
forall a. HasCallStack => [Char] -> a
error
[Char]
"Float indices could not be compared during matrix creation."
where
iI :: Int
iI = forall a b. (RealFrac a, Integral b) => a -> b
round a
i :: Int
jI :: Int
jI = forall a b. (RealFrac a, Integral b) => a -> b
round a
j :: Int
checkEs :: RealFrac a => Int -> [a] -> [a]
checkEs :: forall a. RealFrac a => Int -> [a] -> [a]
checkEs Int
n [a]
es
| forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
es forall a. Eq a => a -> a -> Bool
== Int
nExp = [a]
es
| Bool
otherwise = forall a. HasCallStack => [Char] -> a
error [Char]
eStr
where
nExp :: Int
nExp = forall a b. (RealFrac a, Integral b) => a -> b
round (Int
n Int -> Int -> Double
`choose` Int
2)
eStr :: [Char]
eStr =
[[Char]] -> [Char]
unlines
[ [Char]
"exchFromListlower: the number of exchangeabilities does not match the matrix size",
[Char]
"matrix size: " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> [Char]
show Int
n,
[Char]
"expected number of exchangeabilities: " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> [Char]
show Int
nExp,
[Char]
"received number of exchangeabilities: " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> [Char]
show (forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
es)
]
exchFromListLower :: (RealFrac a, Container Vector a) => Int -> [a] -> Matrix a
exchFromListLower :: forall a.
(RealFrac a, Container Vector a) =>
Int -> [a] -> Matrix a
exchFromListLower Int
n [a]
es = forall d f (c :: * -> *) e. Build d f c e => d -> f -> c e
build (Int
n, Int
n) (forall a. RealFrac a => [a] -> a -> a -> a
fromListBuilderLower (forall a. RealFrac a => Int -> [a] -> [a]
checkEs Int
n [a]
es))
exchFromListUpper :: (RealFrac a, Container Vector a) => Int -> [a] -> Matrix a
exchFromListUpper :: forall a.
(RealFrac a, Container Vector a) =>
Int -> [a] -> Matrix a
exchFromListUpper Int
n [a]
es = forall d f (c :: * -> *) e. Build d f c e => d -> f -> c e
build (Int
n, Int
n) (forall a. RealFrac a => Int -> [a] -> a -> a -> a
fromListBuilderUpper Int
n (forall a. RealFrac a => Int -> [a] -> [a]
checkEs Int
n [a]
es))