-- |
-- Module      :  ELynx.Data.MarkovProcess.CXXModels
-- Description :  C10 to C60 models
-- Copyright   :  (c) Dominik Schrempf 2021
-- License     :  GPL-3.0-or-later
--
-- Maintainer  :  dominik.schrempf@gmail.com
-- Stability   :  unstable
-- Portability :  portable
--
-- Creation date: Tue Feb 26 16:44:33 2019.
--
-- Quang, BL. S., Gascuel, O. & Lartillot, N. Empirical profile mixture models for
-- phylogenetic reconstruction. Bioinformatics 24, 2317–2323 (2008).
--
-- XXX: For now, I only provide Poisson exchangeabilities.
module ELynx.Data.MarkovProcess.CXXModels
  ( cxx,
  )
where

import qualified Data.Vector as V
import ELynx.Data.MarkovProcess.AminoAcid
import ELynx.Data.MarkovProcess.CXXModelsData
import qualified ELynx.Data.MarkovProcess.MixtureModel as M
import ELynx.Data.MarkovProcess.RateMatrix
import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S

-- | Create CXX model with given number of components and probably with custom
-- weights.
cxx :: Int -> Maybe [M.Weight] -> M.MixtureModel
cxx :: Int -> Maybe [Weight] -> MixtureModel
cxx Int
10 (Just [Weight]
ws) = [Weight] -> MixtureModel
c10CustomWeights [Weight]
ws
cxx Int
20 (Just [Weight]
ws) = [Weight] -> MixtureModel
c20CustomWeights [Weight]
ws
cxx Int
30 (Just [Weight]
ws) = [Weight] -> MixtureModel
c30CustomWeights [Weight]
ws
cxx Int
40 (Just [Weight]
ws) = [Weight] -> MixtureModel
c40CustomWeights [Weight]
ws
cxx Int
50 (Just [Weight]
ws) = [Weight] -> MixtureModel
c50CustomWeights [Weight]
ws
cxx Int
60 (Just [Weight]
ws) = [Weight] -> MixtureModel
c60CustomWeights [Weight]
ws
cxx Int
10 Maybe [Weight]
Nothing = MixtureModel
c10
cxx Int
20 Maybe [Weight]
Nothing = MixtureModel
c20
cxx Int
30 Maybe [Weight]
Nothing = MixtureModel
c30
cxx Int
40 Maybe [Weight]
Nothing = MixtureModel
c40
cxx Int
50 Maybe [Weight]
Nothing = MixtureModel
c50
cxx Int
60 Maybe [Weight]
Nothing = MixtureModel
c60
cxx Int
n Maybe [Weight]
_ =
  [Char] -> MixtureModel
forall a. HasCallStack => [Char] -> a
error ([Char] -> MixtureModel) -> [Char] -> MixtureModel
forall a b. (a -> b) -> a -> b
$ [Char]
"cxx: cannot create CXX model with " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show Int
n [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
" components."

-- | C10 model.
c10 :: M.MixtureModel
c10 :: MixtureModel
c10 = [Weight] -> [StationaryDistribution] -> MixtureModel
cxxFromStatDistsAndWeights [Weight]
c10Weights [StationaryDistribution]
c10StatDists

-- | C20 model.
c20 :: M.MixtureModel
c20 :: MixtureModel
c20 = [Weight] -> [StationaryDistribution] -> MixtureModel
cxxFromStatDistsAndWeights [Weight]
c20Weights [StationaryDistribution]
c20StatDists

-- | C30 model.
c30 :: M.MixtureModel
c30 :: MixtureModel
c30 = [Weight] -> [StationaryDistribution] -> MixtureModel
cxxFromStatDistsAndWeights [Weight]
c30Weights [StationaryDistribution]
c30StatDists

-- | C40 model.
c40 :: M.MixtureModel
c40 :: MixtureModel
c40 = [Weight] -> [StationaryDistribution] -> MixtureModel
cxxFromStatDistsAndWeights [Weight]
c40Weights [StationaryDistribution]
c40StatDists

-- | C50 model.
c50 :: M.MixtureModel
c50 :: MixtureModel
c50 = [Weight] -> [StationaryDistribution] -> MixtureModel
cxxFromStatDistsAndWeights [Weight]
c50Weights [StationaryDistribution]
c50StatDists

-- | C60 model.
c60 :: M.MixtureModel
c60 :: MixtureModel
c60 = [Weight] -> [StationaryDistribution] -> MixtureModel
cxxFromStatDistsAndWeights [Weight]
c60Weights [StationaryDistribution]
c60StatDists

-- | C10 model with custom weights.
c10CustomWeights :: [M.Weight] -> M.MixtureModel
c10CustomWeights :: [Weight] -> MixtureModel
c10CustomWeights [Weight]
ws
  | [Weight] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Weight]
ws Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
10 = [Weight] -> [StationaryDistribution] -> MixtureModel
cxxFromStatDistsAndWeights [Weight]
ws [StationaryDistribution]
c10StatDists
  | Bool
otherwise = [Char] -> MixtureModel
forall a. HasCallStack => [Char] -> a
error [Char]
"Number of weights does not match C10 model."

-- | C20 model with custom weights.
c20CustomWeights :: [M.Weight] -> M.MixtureModel
c20CustomWeights :: [Weight] -> MixtureModel
c20CustomWeights [Weight]
ws
  | [Weight] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Weight]
ws Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
20 = [Weight] -> [StationaryDistribution] -> MixtureModel
cxxFromStatDistsAndWeights [Weight]
ws [StationaryDistribution]
c20StatDists
  | Bool
otherwise = [Char] -> MixtureModel
forall a. HasCallStack => [Char] -> a
error [Char]
"Number of weights does not match C20 model."

-- | C30 model with custom weights.
c30CustomWeights :: [M.Weight] -> M.MixtureModel
c30CustomWeights :: [Weight] -> MixtureModel
c30CustomWeights [Weight]
ws
  | [Weight] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Weight]
ws Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
30 = [Weight] -> [StationaryDistribution] -> MixtureModel
cxxFromStatDistsAndWeights [Weight]
ws [StationaryDistribution]
c30StatDists
  | Bool
otherwise = [Char] -> MixtureModel
forall a. HasCallStack => [Char] -> a
error [Char]
"Number of weights does not match C30 model."

-- | C40 model with custom weights.
c40CustomWeights :: [M.Weight] -> M.MixtureModel
c40CustomWeights :: [Weight] -> MixtureModel
c40CustomWeights [Weight]
ws
  | [Weight] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Weight]
ws Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
40 = [Weight] -> [StationaryDistribution] -> MixtureModel
cxxFromStatDistsAndWeights [Weight]
ws [StationaryDistribution]
c40StatDists
  | Bool
otherwise = [Char] -> MixtureModel
forall a. HasCallStack => [Char] -> a
error [Char]
"Number of weights does not match C40 model."

-- | C50 model with custom weights.
c50CustomWeights :: [M.Weight] -> M.MixtureModel
c50CustomWeights :: [Weight] -> MixtureModel
c50CustomWeights [Weight]
ws
  | [Weight] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Weight]
ws Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
50 = [Weight] -> [StationaryDistribution] -> MixtureModel
cxxFromStatDistsAndWeights [Weight]
ws [StationaryDistribution]
c50StatDists
  | Bool
otherwise = [Char] -> MixtureModel
forall a. HasCallStack => [Char] -> a
error [Char]
"Number of weights does not match C50 model."

-- | C60 model with custom weights.
c60CustomWeights :: [M.Weight] -> M.MixtureModel
c60CustomWeights :: [Weight] -> MixtureModel
c60CustomWeights [Weight]
ws
  | [Weight] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Weight]
ws Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
60 = [Weight] -> [StationaryDistribution] -> MixtureModel
cxxFromStatDistsAndWeights [Weight]
ws [StationaryDistribution]
c60StatDists
  | Bool
otherwise = [Char] -> MixtureModel
forall a. HasCallStack => [Char] -> a
error [Char]
"Number of weights does not match C60 model."

cxxName :: Int -> String
cxxName :: Int -> [Char]
cxxName Int
nComps = Char
'C' Char -> [Char] -> [Char]
forall a. a -> [a] -> [a]
: Int -> [Char]
forall a. Show a => a -> [Char]
show Int
nComps

componentName :: Int -> Int -> String
componentName :: Int -> Int -> [Char]
componentName Int
nComps Int
comp = Int -> [Char]
cxxName Int
nComps [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"; component " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show Int
comp

-- Keep in mind, that when using different exchangeabilities, I have to decide
-- about global or local normalization.
cxxSubstitutionModelFromStatDist ::
  Int -> Int -> StationaryDistribution -> S.SubstitutionModel
cxxSubstitutionModelFromStatDist :: Int -> Int -> StationaryDistribution -> SubstitutionModel
cxxSubstitutionModelFromStatDist Int
nComps Int
comp StationaryDistribution
d =
  Maybe [Char] -> StationaryDistribution -> SubstitutionModel
poissonCustom
    ([Char] -> Maybe [Char]
forall a. a -> Maybe a
Just [Char]
name)
    (StationaryDistribution -> StationaryDistribution
normalizeSD StationaryDistribution
d)
  where
    name :: [Char]
name = Int -> Int -> [Char]
componentName Int
nComps Int
comp

cxxSubstitutionModelsFromStatDists ::
  [StationaryDistribution] -> [S.SubstitutionModel]
cxxSubstitutionModelsFromStatDists :: [StationaryDistribution] -> [SubstitutionModel]
cxxSubstitutionModelsFromStatDists [StationaryDistribution]
ds =
  (Int -> StationaryDistribution -> SubstitutionModel)
-> [Int] -> [StationaryDistribution] -> [SubstitutionModel]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith
    (Int -> Int -> StationaryDistribution -> SubstitutionModel
cxxSubstitutionModelFromStatDist Int
nComp)
    [Int
1 ..]
    [StationaryDistribution]
ds
  where
    nComp :: Int
nComp = [StationaryDistribution] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [StationaryDistribution]
ds

-- XXX: The use of `Data.List.NonEmpty.fromList` is daring, but since this
-- function is not exported and only applied to predefined non-empty lists, it
-- should be OK.
cxxFromStatDistsAndWeights ::
  [M.Weight] -> [StationaryDistribution] -> M.MixtureModel
cxxFromStatDistsAndWeights :: [Weight] -> [StationaryDistribution] -> MixtureModel
cxxFromStatDistsAndWeights [Weight]
ws [StationaryDistribution]
ds =
  [Char] -> Vector Weight -> Vector SubstitutionModel -> MixtureModel
M.fromSubstitutionModels
    (Int -> [Char]
cxxName Int
n)
    ([Weight] -> Vector Weight
forall a. [a] -> Vector a
V.fromList [Weight]
ws)
    Vector SubstitutionModel
sms
  where
    n :: Int
n = [StationaryDistribution] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [StationaryDistribution]
ds
    sms :: Vector SubstitutionModel
sms = [SubstitutionModel] -> Vector SubstitutionModel
forall a. [a] -> Vector a
V.fromList ([SubstitutionModel] -> Vector SubstitutionModel)
-> [SubstitutionModel] -> Vector SubstitutionModel
forall a b. (a -> b) -> a -> b
$ [StationaryDistribution] -> [SubstitutionModel]
cxxSubstitutionModelsFromStatDists [StationaryDistribution]
ds