{-# language ViewPatterns, ScopedTypeVariables, MultiWayIf, FlexibleContexts #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Algorithm.SRTree.ConfidenceIntervals 
-- Copyright   :  (c) Fabricio Olivetti 2021 - 2024
-- License     :  BSD3
-- Maintainer  :  fabricio.olivetti@gmail.com
-- Stability   :  experimental
-- Portability :  ConstraintKinds
--
-- Functions to optimize the parameters of an expression.
--
-----------------------------------------------------------------------------
module Algorithm.SRTree.ConfidenceIntervals where

import qualified Data.Massiv.Array as A
import Data.Massiv.Array (Ix2(..), (*.), (!+!), (!*!))
import Data.Massiv.Array.Numeric ( identityMatrix )
import Statistics.Distribution ( ContDistr(quantile) )
import Statistics.Distribution.StudentT ( studentT )
import Statistics.Distribution.FDistribution ( fDistribution )
import qualified Data.Vector.Storable as VS
import Data.SRTree
import Data.SRTree.Eval
import Data.SRTree.Recursion ( cata )
import Algorithm.SRTree.Likelihoods
import Algorithm.SRTree.Opt
    ( minimizeNLLNonUnique, minimizeNLLWithFixedParam )
import Data.List ( sortOn, nubBy )
import Data.Maybe ( fromMaybe )
import Algorithm.SRTree.NonlinearOpt
import Algorithm.Massiv.Utils
import System.IO.Unsafe ( unsafePerformIO )
import Control.Monad.Catch ( catch )

import Debug.Trace ( trace, traceShow )

-- | profile likelihood algorithms: Bates (classical), ODE (faster), Constrained (fastest)
-- The Constrained approach returns only the endpoints.
data PType = Bates | ODE | Constrained deriving (Ix1 -> PType -> ShowS
[PType] -> ShowS
PType -> String
(Ix1 -> PType -> ShowS)
-> (PType -> String) -> ([PType] -> ShowS) -> Show PType
forall a.
(Ix1 -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Ix1 -> PType -> ShowS
showsPrec :: Ix1 -> PType -> ShowS
$cshow :: PType -> String
show :: PType -> String
$cshowList :: [PType] -> ShowS
showList :: [PType] -> ShowS
Show, ReadPrec [PType]
ReadPrec PType
Ix1 -> ReadS PType
ReadS [PType]
(Ix1 -> ReadS PType)
-> ReadS [PType]
-> ReadPrec PType
-> ReadPrec [PType]
-> Read PType
forall a.
(Ix1 -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
$creadsPrec :: Ix1 -> ReadS PType
readsPrec :: Ix1 -> ReadS PType
$creadList :: ReadS [PType]
readList :: ReadS [PType]
$creadPrec :: ReadPrec PType
readPrec :: ReadPrec PType
$creadListPrec :: ReadPrec [PType]
readListPrec :: ReadPrec [PType]
Read)

-- | Confidence Interval using Laplace approximation or profile likelihood.
data CIType = Laplace BasicStats | Profile BasicStats [ProfileT]

-- | Basic stats of the data: covariance of parameters, correlation, standard errors 
data BasicStats = MkStats { BasicStats -> SRMatrix
_cov    :: SRMatrix
                          , BasicStats -> SRMatrix
_corr   :: SRMatrix
                          , BasicStats -> PVector
_stdErr :: PVector
                          } deriving (BasicStats -> BasicStats -> Bool
(BasicStats -> BasicStats -> Bool)
-> (BasicStats -> BasicStats -> Bool) -> Eq BasicStats
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: BasicStats -> BasicStats -> Bool
== :: BasicStats -> BasicStats -> Bool
$c/= :: BasicStats -> BasicStats -> Bool
/= :: BasicStats -> BasicStats -> Bool
Eq, Ix1 -> BasicStats -> ShowS
[BasicStats] -> ShowS
BasicStats -> String
(Ix1 -> BasicStats -> ShowS)
-> (BasicStats -> String)
-> ([BasicStats] -> ShowS)
-> Show BasicStats
forall a.
(Ix1 -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Ix1 -> BasicStats -> ShowS
showsPrec :: Ix1 -> BasicStats -> ShowS
$cshow :: BasicStats -> String
show :: BasicStats -> String
$cshowList :: [BasicStats] -> ShowS
showList :: [BasicStats] -> ShowS
Show)

-- | a confience interval is composed of the point estimate (`est_`), lower bound (`_lower_`)
-- and upper bound (`upper_`)
data CI = CI { CI -> Double
est_   :: Double
             , CI -> Double
lower_ :: Double
             , CI -> Double
upper_ :: Double
             } deriving (CI -> CI -> Bool
(CI -> CI -> Bool) -> (CI -> CI -> Bool) -> Eq CI
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: CI -> CI -> Bool
== :: CI -> CI -> Bool
$c/= :: CI -> CI -> Bool
/= :: CI -> CI -> Bool
Eq, Ix1 -> CI -> ShowS
[CI] -> ShowS
CI -> String
(Ix1 -> CI -> ShowS)
-> (CI -> String) -> ([CI] -> ShowS) -> Show CI
forall a.
(Ix1 -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Ix1 -> CI -> ShowS
showsPrec :: Ix1 -> CI -> ShowS
$cshow :: CI -> String
show :: CI -> String
$cshowList :: [CI] -> ShowS
showList :: [CI] -> ShowS
Show, ReadPrec [CI]
ReadPrec CI
Ix1 -> ReadS CI
ReadS [CI]
(Ix1 -> ReadS CI)
-> ReadS [CI] -> ReadPrec CI -> ReadPrec [CI] -> Read CI
forall a.
(Ix1 -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
$creadsPrec :: Ix1 -> ReadS CI
readsPrec :: Ix1 -> ReadS CI
$creadList :: ReadS [CI]
readList :: ReadS [CI]
$creadPrec :: ReadPrec CI
readPrec :: ReadPrec CI
$creadListPrec :: ReadPrec [CI]
readListPrec :: ReadPrec [CI]
Read)

--  | A profile likelihood is composed of a vector of tau values that traces the likelihood, 
--  the matrix of thetas for each profile, the local optima, and two splines that converts 
--  taus to theta and vice-versa. 
data ProfileT = ProfileT { ProfileT -> PVector
_taus      :: PVector
                         , ProfileT -> SRMatrix
_thetas    :: SRMatrix
                         , ProfileT -> Double
_opt       :: Double
                         , ProfileT -> Double -> Double
_tau2theta :: Double -> Double
                         , ProfileT -> Double -> Double
_theta2tau :: Double -> Double
                         }

-- shows the CI with n places 
showCI :: Int -> CI -> String
showCI :: Ix1 -> CI -> String
showCI Ix1
n (CI Double
x Double
l Double
h) = Double -> String
forall a. Show a => a -> String
show (Double -> Double
rnd Double
l) String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String
" <= " String -> ShowS
forall a. Semigroup a => a -> a -> a
<> Double -> String
forall a. Show a => a -> String
show (Double -> Double
rnd Double
x) String -> ShowS
forall a. Semigroup a => a -> a -> a
<> String
" <= " String -> ShowS
forall a. Semigroup a => a -> a -> a
<> Double -> String
forall a. Show a => a -> String
show (Double -> Double
rnd Double
h)
  where
      rnd :: Double -> Double
rnd = (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
10Double -> Ix1 -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Ix1
n) (Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer -> Double) -> (Double -> Integer) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round) (Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
10Double -> Ix1 -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Ix1
n)
printCI :: Int -> CI -> IO ()
printCI :: Ix1 -> CI -> IO ()
printCI Ix1
n = String -> IO ()
putStrLn (String -> IO ()) -> (CI -> String) -> CI -> IO ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Ix1 -> CI -> String
showCI Ix1
n

-- | Calculates the confidence interval of the parameters using 
-- Laplace approximation or Profile likelihood
paramCI :: CIType -> Int -> PVector -> Double -> [CI]
paramCI :: CIType -> Ix1 -> PVector -> Double -> [CI]
paramCI (Laplace BasicStats
stats) Ix1
nSamples PVector
theta Double
alpha = (Double -> Double -> Double -> CI)
-> [Double] -> [Double] -> [Double] -> [CI]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 Double -> Double -> Double -> CI
CI (PVector -> [Double]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList PVector
theta) [Double]
lows [Double]
highs
  where
    -- the Laplace approximation is theta +/- t(1-alpha/2) * standard error 
    (A.Sz Ix1
k) = PVector -> Sz Ix1
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
A.size PVector
theta
    t :: Double
t        = StudentT -> Double -> Double
forall d. ContDistr d => d -> Double -> Double
quantile (Double -> StudentT
studentT (Double -> StudentT) -> (Ix1 -> Double) -> Ix1 -> StudentT
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Ix1 -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Ix1 -> StudentT) -> Ix1 -> StudentT
forall a b. (a -> b) -> a -> b
$ Ix1
nSamples Ix1 -> Ix1 -> Ix1
forall a. Num a => a -> a -> a
- Ix1
k) (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2.0)
    stdErr :: PVector
stdErr   = BasicStats -> PVector
_stdErr BasicStats
stats
    lows :: [Double]
lows     = Array D Ix1 Double -> [Double]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList (Array D Ix1 Double -> [Double]) -> Array D Ix1 Double -> [Double]
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double)
-> PVector -> Array D Ix1 Double -> Array D Ix1 Double
forall ix r1 e1 r2 e2 e.
(Index ix, Source r1 e1, Source r2 e2) =>
(e1 -> e2 -> e) -> Array r1 ix e1 -> Array r2 ix e2 -> Array D ix e
A.zipWith (-) PVector
theta (Array D Ix1 Double -> Array D Ix1 Double)
-> Array D Ix1 Double -> Array D Ix1 Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double) -> PVector -> Array D Ix1 Double
forall ix r e' e.
(Index ix, Source r e') =>
(e' -> e) -> Array r ix e' -> Array D ix e
A.map (Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t) PVector
stdErr
    highs :: [Double]
highs    = Array D Ix1 Double -> [Double]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList (Array D Ix1 Double -> [Double]) -> Array D Ix1 Double -> [Double]
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double)
-> PVector -> Array D Ix1 Double -> Array D Ix1 Double
forall ix r1 e1 r2 e2 e.
(Index ix, Source r1 e1, Source r2 e2) =>
(e1 -> e2 -> e) -> Array r1 ix e1 -> Array r2 ix e2 -> Array D ix e
A.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) PVector
theta (Array D Ix1 Double -> Array D Ix1 Double)
-> Array D Ix1 Double -> Array D Ix1 Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double) -> PVector -> Array D Ix1 Double
forall ix r e' e.
(Index ix, Source r e') =>
(e' -> e) -> Array r ix e' -> Array D ix e
A.map (Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t) PVector
stdErr

paramCI (Profile BasicStats
stats [ProfileT]
profiles) Ix1
nSamples PVector
_ Double
alpha = (Double -> Double -> Double -> CI)
-> [Double] -> [Double] -> [Double] -> [CI]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 Double -> Double -> Double -> CI
CI [Double]
theta [Double]
lows [Double]
highs
  where
    -- for the profile likelihood we use the square root of the F-distribution with (1-alpha)
    k :: Ix1
k        = [Double] -> Ix1
forall a. [a] -> Ix1
forall (t :: * -> *) a. Foldable t => t a -> Ix1
length [Double]
theta
    t :: Double
t        = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ FDistribution -> Double -> Double
forall d. ContDistr d => d -> Double -> Double
quantile (Ix1 -> Ix1 -> FDistribution
fDistribution Ix1
k (Ix1 -> Ix1
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Ix1 -> Ix1) -> Ix1 -> Ix1
forall a b. (a -> b) -> a -> b
$ Ix1
nSamples Ix1 -> Ix1 -> Ix1
forall a. Num a => a -> a -> a
- Ix1
k)) (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha)
    stdErr :: PVector
stdErr   = BasicStats -> PVector
_stdErr BasicStats
stats
    lows :: [Double]
lows     = (ProfileT -> Double) -> [ProfileT] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (ProfileT -> Double -> Double
`_tau2theta` (-Double
t)) [ProfileT]
profiles
    highs :: [Double]
highs    = (ProfileT -> Double) -> [ProfileT] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (ProfileT -> Double -> Double
`_tau2theta` Double
t) [ProfileT]
profiles
    theta :: [Double]
theta    = (ProfileT -> Double) -> [ProfileT] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ProfileT -> Double
_opt [ProfileT]
profiles

-- | calculates the prediction confidence interval using Laplace approximation or profile likelihood. 
--
predictionCI :: CIType -> Distribution -> (SRMatrix -> PVector) -> (SRMatrix -> [PVector]) -> (CI -> PVector -> Fix SRTree -> (Double -> Double, Double)) -> SRMatrix -> Fix SRTree -> PVector -> Double -> [CI] -> [CI]
predictionCI :: CIType
-> Distribution
-> (SRMatrix -> PVector)
-> (SRMatrix -> [PVector])
-> (CI -> PVector -> Fix SRTree -> (Double -> Double, Double))
-> SRMatrix
-> Fix SRTree
-> PVector
-> Double
-> [CI]
-> [CI]
predictionCI (Laplace BasicStats
stats) Distribution
_ SRMatrix -> PVector
predFun SRMatrix -> [PVector]
jacFun CI -> PVector -> Fix SRTree -> (Double -> Double, Double)
_ SRMatrix
xss Fix SRTree
tree PVector
theta Double
alpha [CI]
_ = (Double -> Double -> Double -> CI)
-> [Double] -> [Double] -> [Double] -> [CI]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 Double -> Double -> Double -> CI
CI [Double]
yhat [Double]
lows [Double]
highs
  where
    yhat :: [Double]
yhat     = PVector -> [Double]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList (PVector -> [Double]) -> PVector -> [Double]
forall a b. (a -> b) -> a -> b
$ SRMatrix -> PVector
predFun SRMatrix
xss
    jac' :: A.Matrix A.S Double
    jac' :: SRMatrix
jac'     = Comp -> [ListItem Ix2 Double] -> SRMatrix
forall r ix e.
(HasCallStack, Ragged L ix e, Manifest r e) =>
Comp -> [ListItem ix e] -> Array r ix e
A.fromLists' Comp
compMode ([ListItem Ix2 Double] -> SRMatrix)
-> [ListItem Ix2 Double] -> SRMatrix
forall a b. (a -> b) -> a -> b
$ (PVector -> ListItem Ix2 Double)
-> [PVector] -> [ListItem Ix2 Double]
forall a b. (a -> b) -> [a] -> [b]
map PVector -> [Double]
PVector -> ListItem Ix2 Double
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList ([PVector] -> [ListItem Ix2 Double])
-> [PVector] -> [ListItem Ix2 Double]
forall a b. (a -> b) -> a -> b
$ SRMatrix -> [PVector]
jacFun SRMatrix
xss
    jac :: [PVector]
    jac :: [PVector]
jac      = Array D Ix1 PVector -> [PVector]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList (Array D Ix1 PVector -> [PVector])
-> Array D Ix1 PVector -> [PVector]
forall a b. (a -> b) -> a -> b
$ SRMatrix -> Array D Ix1 (Array S (Lower Ix2) Double)
forall r ix e.
(Index ix, Index (Lower ix), Source r e) =>
Array r ix e -> Array D Ix1 (Array r (Lower ix) e)
A.outerSlices (SRMatrix -> Array D Ix1 (Array S (Lower Ix2) Double))
-> SRMatrix -> Array D Ix1 (Array S (Lower Ix2) Double)
forall a b. (a -> b) -> a -> b
$ S -> Array D Ix2 Double -> SRMatrix
forall r e r' ix.
(Manifest r e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
A.computeAs S
A.S (Array D Ix2 Double -> SRMatrix) -> Array D Ix2 Double -> SRMatrix
forall a b. (a -> b) -> a -> b
$ SRMatrix -> Array D Ix2 Double
forall r e. Source r e => Matrix r e -> Matrix D e
A.transpose SRMatrix
jac'
    n :: Ix1
n        = [Double] -> Ix1
forall a. [a] -> Ix1
forall (t :: * -> *) a. Foldable t => t a -> Ix1
length [Double]
yhat
    (A.Sz Ix1
k) = PVector -> Sz Ix1
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
A.size PVector
theta
    t :: Double
t        = StudentT -> Double -> Double
forall d. ContDistr d => d -> Double -> Double
quantile (Double -> StudentT
studentT (Double -> StudentT) -> (Ix1 -> Double) -> Ix1 -> StudentT
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Ix1 -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Ix1 -> StudentT) -> Ix1 -> StudentT
forall a b. (a -> b) -> a -> b
$ Ix1
n Ix1 -> Ix1 -> Ix1
forall a. Num a => a -> a -> a
- Ix1
k) (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2.0)
    covs :: [PVector]
covs     = Array D Ix1 PVector -> [PVector]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList (Array D Ix1 PVector -> [PVector])
-> Array D Ix1 PVector -> [PVector]
forall a b. (a -> b) -> a -> b
$ SRMatrix -> Array D Ix1 (Array S (Lower Ix2) Double)
forall r ix e.
(Index ix, Index (Lower ix), Source r e) =>
Array r ix e -> Array D Ix1 (Array r (Lower ix) e)
A.outerSlices (SRMatrix -> Array D Ix1 (Array S (Lower Ix2) Double))
-> SRMatrix -> Array D Ix1 (Array S (Lower Ix2) Double)
forall a b. (a -> b) -> a -> b
$ BasicStats -> SRMatrix
_cov BasicStats
stats
    lows :: [Double]
lows     = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
yhat ([Double] -> [Double]) -> [Double] -> [Double]
forall a b. (a -> b) -> a -> b
$ (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t) [Double]
resStdErr
    highs :: [Double]
highs    = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) [Double]
yhat ([Double] -> [Double]) -> [Double] -> [Double]
forall a b. (a -> b) -> a -> b
$ (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t) [Double]
resStdErr

    getResStdError :: PVector -> Double
getResStdError PVector
row = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ PVector -> PVector -> Double
forall r e.
(Numeric r e, Source r e) =>
Vector r e -> Vector r e -> e
(A.!.!) PVector
row (PVector -> Double) -> PVector -> Double
forall a b. (a -> b) -> a -> b
$ Comp -> [Double] -> PVector
forall r e. Manifest r e => Comp -> [e] -> Vector r e
A.fromList Comp
compMode ([Double] -> PVector) -> [Double] -> PVector
forall a b. (a -> b) -> a -> b
$ (PVector -> Double) -> [PVector] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (PVector
row PVector -> PVector -> Double
forall r e.
(Numeric r e, Source r e) =>
Vector r e -> Vector r e -> e
A.!.!) [PVector]
covs
    resStdErr :: [Double]
resStdErr          = (PVector -> Double) -> [PVector] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map PVector -> Double
getResStdError [PVector]
jac

predictionCI (Profile BasicStats
_ [ProfileT]
_) Distribution
dist SRMatrix -> PVector
predFun SRMatrix -> [PVector]
_ CI -> PVector -> Fix SRTree -> (Double -> Double, Double)
profFun SRMatrix
xss Fix SRTree
tree PVector
theta Double
alpha [CI]
estPIs = (CI -> Double -> PVector -> CI)
-> [CI] -> [Double] -> [PVector] -> [CI]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 CI -> Double -> PVector -> CI
f [CI]
estPIs [Double]
yhat ([PVector] -> [CI]) -> [PVector] -> [CI]
forall a b. (a -> b) -> a -> b
$ Ix1 -> [PVector] -> [PVector]
forall a. Ix1 -> [a] -> [a]
take Ix1
10 [PVector]
xss'
  where
    yhat :: [Double]
yhat     = PVector -> [Double]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList (PVector -> [Double]) -> PVector -> [Double]
forall a b. (a -> b) -> a -> b
$ SRMatrix -> PVector
predFun SRMatrix
xss
    theta' :: Vector Double
theta'   = PVector -> Vector Double
forall ix e. Index ix => Array S ix e -> Vector e
A.toStorableVector PVector
theta

    t :: Double
t        = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ FDistribution -> Double -> Double
forall d. ContDistr d => d -> Double -> Double
quantile (Ix1 -> Ix1 -> FDistribution
fDistribution Ix1
k (Ix1 -> Ix1
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Ix1 -> Ix1) -> Ix1 -> Ix1
forall a b. (a -> b) -> a -> b
$ Ix1
n Ix1 -> Ix1 -> Ix1
forall a. Num a => a -> a -> a
- Ix1
k)) (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha)
    (A.Sz Ix1
k) = PVector -> Sz Ix1
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
A.size PVector
theta
    n :: Ix1
n        = [Double] -> Ix1
forall a. [a] -> Ix1
forall (t :: * -> *) a. Foldable t => t a -> Ix1
length [Double]
yhat

    theta0 :: Fix SRTree
theta0  = Distribution -> Fix SRTree -> Fix SRTree
calcTheta0 Distribution
dist Fix SRTree
tree
    xss' :: [PVector]
xss'    = Array D Ix1 PVector -> [PVector]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList (Array D Ix1 PVector -> [PVector])
-> Array D Ix1 PVector -> [PVector]
forall a b. (a -> b) -> a -> b
$ SRMatrix -> Array D Ix1 (Array S (Lower Ix2) Double)
forall r ix e.
(Index ix, Index (Lower ix), Source r e) =>
Array r ix e -> Array D Ix1 (Array r (Lower ix) e)
A.outerSlices SRMatrix
xss

    f :: CI -> Double -> PVector -> CI
f CI
estPI Double
yh PVector
xs =
              let t' :: Fix SRTree
t'            = Fix SRTree -> Fix SRTree -> Fix SRTree
replaceParam0 Fix SRTree
tree (Fix SRTree -> Fix SRTree) -> Fix SRTree -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ PVector -> Fix SRTree -> Fix SRTree
evalVar PVector
xs Fix SRTree
theta0
                  (Double -> Double
spline, Double
yh') = CI -> PVector -> Fix SRTree -> (Double -> Double, Double)
profFun CI
estPI (Comp -> Vector Double -> PVector
forall e. Comp -> Vector e -> Vector S e
A.fromStorableVector Comp
compMode (Vector Double
theta' Vector Double -> [(Ix1, Double)] -> Vector Double
forall a. Storable a => Vector a -> [(Ix1, a)] -> Vector a
VS.// [(Ix1
0, Double
yh)])) Fix SRTree
t'
              in Double -> Double -> Double -> CI
CI Double
yh' (Double -> Double
spline (-Double
t)) (Double -> Double
spline Double
t)

-- inverse function of the distributions 
inverseDist :: Floating p => Distribution -> p -> p
inverseDist :: forall p. Floating p => Distribution -> p -> p
inverseDist Distribution
Gaussian p
y  = p
y
inverseDist Distribution
Bernoulli p
y = p -> p
forall a. Floating a => a -> a
log (p
yp -> p -> p
forall a. Fractional a => a -> a -> a
/(p
1p -> p -> p
forall a. Num a => a -> a -> a
-p
y))
inverseDist Distribution
Poisson p
y   = p -> p
forall a. Floating a => a -> a
log p
y

-- rewrite the tree by fixing theta 0 to optimal value 
replaceParam0 :: Fix SRTree -> Fix SRTree -> Fix SRTree
replaceParam0 :: Fix SRTree -> Fix SRTree -> Fix SRTree
replaceParam0 Fix SRTree
tree Fix SRTree
t0 = (SRTree (Fix SRTree) -> Fix SRTree) -> Fix SRTree -> Fix SRTree
forall (f :: * -> *) a. Functor f => (f a -> a) -> Fix f -> a
cata SRTree (Fix SRTree) -> Fix SRTree
alg Fix SRTree
tree
  where
    alg :: SRTree (Fix SRTree) -> Fix SRTree
alg (Var Ix1
ix)     = SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Ix1 -> SRTree (Fix SRTree)
forall val. Ix1 -> SRTree val
Var Ix1
ix
    alg (Param Ix1
0)    = Fix SRTree
t0
    alg (Param Ix1
ix)   = SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Ix1 -> SRTree (Fix SRTree)
forall val. Ix1 -> SRTree val
Param Ix1
ix
    alg (Const Double
c)    = SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Double -> SRTree (Fix SRTree)
forall val. Double -> SRTree val
Const Double
c
    alg (Uni Function
g Fix SRTree
t)    = SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Function -> Fix SRTree -> SRTree (Fix SRTree)
forall val. Function -> val -> SRTree val
Uni Function
g Fix SRTree
t
    alg (Bin Op
op Fix SRTree
l Fix SRTree
r) = SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Op -> Fix SRTree -> Fix SRTree -> SRTree (Fix SRTree)
forall val. Op -> val -> val -> SRTree val
Bin Op
op Fix SRTree
l Fix SRTree
r

evalVar :: PVector -> Fix SRTree -> Fix SRTree
evalVar :: PVector -> Fix SRTree -> Fix SRTree
evalVar PVector
xs = (SRTree (Fix SRTree) -> Fix SRTree) -> Fix SRTree -> Fix SRTree
forall (f :: * -> *) a. Functor f => (f a -> a) -> Fix f -> a
cata SRTree (Fix SRTree) -> Fix SRTree
alg
  where
    alg :: SRTree (Fix SRTree) -> Fix SRTree
alg (Var Ix1
ix)     = SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Double -> SRTree (Fix SRTree)
forall val. Double -> SRTree val
Const (PVector
xs PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix)
    alg (Param Ix1
ix)   = SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Ix1 -> SRTree (Fix SRTree)
forall val. Ix1 -> SRTree val
Param Ix1
ix
    alg (Const Double
c)    = SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Double -> SRTree (Fix SRTree)
forall val. Double -> SRTree val
Const Double
c
    alg (Uni Function
g Fix SRTree
t)    = SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Function -> Fix SRTree -> SRTree (Fix SRTree)
forall val. Function -> val -> SRTree val
Uni Function
g Fix SRTree
t
    alg (Bin Op
op Fix SRTree
l Fix SRTree
r) = SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Op -> Fix SRTree -> Fix SRTree -> SRTree (Fix SRTree)
forall val. Op -> val -> val -> SRTree val
Bin Op
op Fix SRTree
l Fix SRTree
r

calcTheta0 :: Distribution -> Fix SRTree -> Fix SRTree
calcTheta0 :: Distribution -> Fix SRTree -> Fix SRTree
calcTheta0 Distribution
dist Fix SRTree
tree = case (SRTree (Either (Fix SRTree -> Fix SRTree) (Fix SRTree))
 -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree))
-> Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall (f :: * -> *) a. Functor f => (f a -> a) -> Fix f -> a
cata SRTree (Either (Fix SRTree -> Fix SRTree) (Fix SRTree))
-> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
alg Fix SRTree
tree of
                         Left Fix SRTree -> Fix SRTree
g -> Fix SRTree -> Fix SRTree
g (Fix SRTree -> Fix SRTree) -> Fix SRTree -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Distribution -> Fix SRTree -> Fix SRTree
forall p. Floating p => Distribution -> p -> p
inverseDist Distribution
dist (SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Ix1 -> SRTree (Fix SRTree)
forall val. Ix1 -> SRTree val
Param Ix1
0)
                         Right Fix SRTree
_ -> String -> Fix SRTree
forall a. HasCallStack => String -> a
error String
"No theta0?"
  where
    alg :: SRTree (Either (Fix SRTree -> Fix SRTree) (Fix SRTree))
-> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
alg (Var Ix1
ix)     = Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. b -> Either a b
Right (Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree))
-> Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. (a -> b) -> a -> b
$ SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Ix1 -> SRTree (Fix SRTree)
forall val. Ix1 -> SRTree val
Var Ix1
ix
    alg (Param Ix1
0)    = (Fix SRTree -> Fix SRTree)
-> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. a -> Either a b
Left Fix SRTree -> Fix SRTree
forall a. a -> a
id
    alg (Param Ix1
ix)   = Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. b -> Either a b
Right (Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree))
-> Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. (a -> b) -> a -> b
$ SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Ix1 -> SRTree (Fix SRTree)
forall val. Ix1 -> SRTree val
Param Ix1
ix
    alg (Const Double
c)    = Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. b -> Either a b
Right (Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree))
-> Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. (a -> b) -> a -> b
$ SRTree (Fix SRTree) -> Fix SRTree
forall (f :: * -> *). f (Fix f) -> Fix f
Fix (SRTree (Fix SRTree) -> Fix SRTree)
-> SRTree (Fix SRTree) -> Fix SRTree
forall a b. (a -> b) -> a -> b
$ Double -> SRTree (Fix SRTree)
forall val. Double -> SRTree val
Const Double
c
    alg (Uni Function
g Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
t)    = case Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
t of
                         Left Fix SRTree -> Fix SRTree
f  -> (Fix SRTree -> Fix SRTree)
-> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. a -> Either a b
Left ((Fix SRTree -> Fix SRTree)
 -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree))
-> (Fix SRTree -> Fix SRTree)
-> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. (a -> b) -> a -> b
$ Fix SRTree -> Fix SRTree
f (Fix SRTree -> Fix SRTree)
-> (Fix SRTree -> Fix SRTree) -> Fix SRTree -> Fix SRTree
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Function -> Fix SRTree -> Fix SRTree
forall a. Floating a => Function -> a -> a
evalInverse Function
g
                         Right Fix SRTree
v -> Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. b -> Either a b
Right (Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree))
-> Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. (a -> b) -> a -> b
$ Function -> Fix SRTree -> Fix SRTree
forall a. Floating a => Function -> a -> a
evalFun Function
g Fix SRTree
v
    alg (Bin Op
op Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
l Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
r) = case Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
l of
                         Left Fix SRTree -> Fix SRTree
f   -> case Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
r of
                                       Left  Fix SRTree -> Fix SRTree
_ -> String -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a. HasCallStack => String -> a
error String
"This shouldn't happen!"
                                       Right Fix SRTree
v -> (Fix SRTree -> Fix SRTree)
-> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. a -> Either a b
Left ((Fix SRTree -> Fix SRTree)
 -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree))
-> (Fix SRTree -> Fix SRTree)
-> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. (a -> b) -> a -> b
$ Fix SRTree -> Fix SRTree
f (Fix SRTree -> Fix SRTree)
-> (Fix SRTree -> Fix SRTree) -> Fix SRTree -> Fix SRTree
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Op -> Fix SRTree -> Fix SRTree -> Fix SRTree
forall a. Floating a => Op -> a -> a -> a
invright Op
op Fix SRTree
v
                         Right Fix SRTree
vl -> case Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
r of
                                       Left  Fix SRTree -> Fix SRTree
g -> (Fix SRTree -> Fix SRTree)
-> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. a -> Either a b
Left ((Fix SRTree -> Fix SRTree)
 -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree))
-> (Fix SRTree -> Fix SRTree)
-> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. (a -> b) -> a -> b
$ Fix SRTree -> Fix SRTree
g (Fix SRTree -> Fix SRTree)
-> (Fix SRTree -> Fix SRTree) -> Fix SRTree -> Fix SRTree
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Op -> Fix SRTree -> Fix SRTree -> Fix SRTree
forall a. Floating a => Op -> a -> a -> a
invleft Op
op Fix SRTree
vl
                                       Right Fix SRTree
vr -> Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. b -> Either a b
Right (Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree))
-> Fix SRTree -> Either (Fix SRTree -> Fix SRTree) (Fix SRTree)
forall a b. (a -> b) -> a -> b
$ Op -> Fix SRTree -> Fix SRTree -> Fix SRTree
forall a. Floating a => Op -> a -> a -> a
evalOp Op
op Fix SRTree
vl Fix SRTree
vr

-- calculate the profile likelihood of every parameter 
getAllProfiles :: PType -> Distribution -> Maybe Double -> SRMatrix -> PVector -> Fix SRTree -> PVector -> PVector -> [CI] -> Double -> [ProfileT]
getAllProfiles :: PType
-> Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> PVector
-> [CI]
-> Double
-> [ProfileT]
getAllProfiles PType
ptype Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta PVector
stdErr [CI]
estCIs Double
alpha = [ProfileT] -> [ProfileT]
forall a. [a] -> [a]
reverse (Ix1 -> [ProfileT] -> [ProfileT]
getAll Ix1
0 [])
  where
    (A.Sz Ix1
k)   = PVector -> Sz Ix1
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
A.size PVector
theta
    (A.Sz Ix1
n)   = PVector -> Sz Ix1
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
A.size PVector
ys
    tau_max :: Double
tau_max    = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ FDistribution -> Double -> Double
forall d. ContDistr d => d -> Double -> Double
quantile (Ix1 -> Ix1 -> FDistribution
fDistribution Ix1
k (Ix1
n Ix1 -> Ix1 -> Ix1
forall a. Num a => a -> a -> a
- Ix1
k)) (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.01)
    tau_max' :: Double
tau_max'   = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ FDistribution -> Double -> Double
forall d. ContDistr d => d -> Double -> Double
quantile (Ix1 -> Ix1 -> FDistribution
fDistribution Ix1
k (Ix1
n Ix1 -> Ix1 -> Ix1
forall a. Num a => a -> a -> a
- Ix1
k)) (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha)

    profFun :: Ix1 -> Either PVector ProfileT
profFun Ix1
ix = case PType
ptype of
                    PType
Bates       -> Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
-> Double
-> Ix1
-> Either PVector ProfileT
getProfile      Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta (PVector
stdErr PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix) Double
tau_max Ix1
ix
                    PType
ODE         -> Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
-> CI
-> Double
-> Ix1
-> Either PVector ProfileT
getProfileODE   Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta (PVector
stdErr PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix) ([CI]
estCIs [CI] -> Ix1 -> CI
forall a. HasCallStack => [a] -> Ix1 -> a
!! Ix1
ix) Double
tau_max Ix1
ix
                    PType
Constrained -> Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
-> Double
-> Ix1
-> Either PVector ProfileT
getProfileCnstr Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta (PVector
stdErr PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix) Double
tau_max' Ix1
ix

    getAll :: Ix1 -> [ProfileT] -> [ProfileT]
getAll Ix1
ix [ProfileT]
acc | Ix1
ix Ix1 -> Ix1 -> Bool
forall a. Eq a => a -> a -> Bool
== Ix1
k   = [ProfileT]
acc
                  | Bool
otherwise = case Ix1 -> Either PVector ProfileT
profFun Ix1
ix of
                                  Left PVector
t  -> PType
-> Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> PVector
-> [CI]
-> Double
-> [ProfileT]
getAllProfiles PType
ptype Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
t PVector
stdErr [CI]
estCIs Double
alpha
                                  Right ProfileT
p -> Ix1 -> [ProfileT] -> [ProfileT]
getAll (Ix1
ix Ix1 -> Ix1 -> Ix1
forall a. Num a => a -> a -> a
+ Ix1
1) (ProfileT
p ProfileT -> [ProfileT] -> [ProfileT]
forall a. a -> [a] -> [a]
: [ProfileT]
acc)

-- calculates the profile likelihood of a single parameter 
getProfile :: Distribution
           -> Maybe Double
           -> SRMatrix
           -> PVector
           -> Fix SRTree
           -> PVector
           -> Double
           -> Double
           -> Int
           -> Either PVector ProfileT
getProfile :: Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
-> Double
-> Ix1
-> Either PVector ProfileT
getProfile Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta Double
stdErr_i Double
tau_max Ix1
ix
  | Double
stdErr_i Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0.0 = ProfileT -> Either PVector ProfileT
forall a. a -> Either PVector a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (ProfileT -> Either PVector ProfileT)
-> ProfileT -> Either PVector ProfileT
forall a b. (a -> b) -> a -> b
$ PVector
-> SRMatrix
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> ProfileT
ProfileT (Comp -> [Double] -> PVector
forall r e. Manifest r e => Comp -> [e] -> Vector r e
A.fromList Comp
compMode [-Double
tau_max, Double
tau_max]) (Comp -> [ListItem Ix2 Double] -> SRMatrix
forall r ix e.
(HasCallStack, Ragged L ix e, Manifest r e) =>
Comp -> [ListItem ix e] -> Array r ix e
A.fromLists' Comp
compMode [[Double]
ListItem Ix2 Double
theta', [Double]
ListItem Ix2 Double
theta']) (PVector
theta PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix) (Double -> Double -> Double
forall a b. a -> b -> a
const (PVector
theta PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix)) (Double -> Double -> Double
forall a b. a -> b -> a
const Double
tau_max)
  | Bool
otherwise =
  do negDelta <- Integer
-> Double
-> Double
-> Double
-> ([Double], [PVector])
-> Either PVector ([Double], [PVector])
forall {t}.
(Eq t, Num t) =>
t
-> Double
-> Double
-> Double
-> ([Double], [PVector])
-> Either PVector ([Double], [PVector])
go Integer
kmax (-Double
stdErr_i Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
8) Double
0 Double
1 ([Double], [PVector])
forall a. Monoid a => a
mempty
     posDelta <- go kmax  (stdErr_i / 8) 0 1 p0
     let (A.fromList compMode -> taus, A.fromLists' compMode. map A.toList -> thetas) = negDelta <> posDelta
         (tau2theta, theta2tau)                       = createSplines taus thetas stdErr_i tau_max ix
     pure $ ProfileT taus thetas optTh tau2theta theta2tau
  where
    theta' :: [Double]
theta'    = PVector -> [Double]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList PVector
theta
    p0 :: ([Double], [PVector])
p0        = ([Double
0], [PVector
theta_opt])
    kmax :: Integer
kmax      = Integer
300
    nll_opt :: Double
nll_opt   = Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
nll Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta_opt
    theta_opt :: PVector
theta_opt = (PVector, Double) -> PVector
forall a b. (a, b) -> a
fst ((PVector, Double) -> PVector) -> (PVector, Double) -> PVector
forall a b. (a -> b) -> a -> b
$ Distribution
-> Maybe Double
-> Ix1
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (PVector, Double)
minimizeNLLNonUnique Distribution
dist Maybe Double
mSErr Ix1
100 SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta
    optTh :: Double
optTh     = PVector
theta_opt PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix
    minimizer :: PVector -> PVector
minimizer = Distribution
-> Maybe Double
-> Ix1
-> SRMatrix
-> PVector
-> Fix SRTree
-> Ix1
-> PVector
-> PVector
minimizeNLLWithFixedParam Distribution
dist Maybe Double
mSErr Ix1
100 SRMatrix
xss PVector
ys Fix SRTree
tree Ix1
ix

    -- after k iterations, interpolates to the endpoint
    go :: t
-> Double
-> Double
-> Double
-> ([Double], [PVector])
-> Either PVector ([Double], [PVector])
go t
0 Double
delta Double
_ Double
_         ([Double], [PVector])
acc = ([Double], [PVector]) -> Either PVector ([Double], [PVector])
forall a b. b -> Either a b
Right ([Double], [PVector])
acc
    go t
k Double
delta Double
t Double
inv_slope acc :: ([Double], [PVector])
acc@([Double]
taus, [PVector]
thetas)
      | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
inv_slope     = ([Double], [PVector]) -> Either PVector ([Double], [PVector])
forall a b. b -> Either a b
Right ([Double], [PVector])
acc    -- stop since we cannot move forward on discontinuity
      | Double
nll_cond Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
nll_opt  = PVector -> Either PVector ([Double], [PVector])
forall a b. a -> Either a b
Left PVector
theta_t -- found a better optima
      | Double -> Double
forall a. Num a => a -> a
abs Double
tau Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
tau_max   = ([Double], [PVector]) -> Either PVector ([Double], [PVector])
forall a b. b -> Either a b
Right ([Double], [PVector])
acc'   -- we reached the endpoint
      | Bool
otherwise           = t
-> Double
-> Double
-> Double
-> ([Double], [PVector])
-> Either PVector ([Double], [PVector])
go (t
kt -> t -> t
forall a. Num a => a -> a -> a
-t
1) Double
delta (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
inv_slope) Double
inv_slope' ([Double], [PVector])
acc'
      where
        t_delta :: Double
t_delta     = (PVector
theta_opt PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
delta Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
inv_slope)
        theta_delta :: PVector
theta_delta = PVector -> [(Ix1, Double)] -> PVector
updateS PVector
theta_opt [(Ix1
ix, Double
t_delta)]
        theta_t :: PVector
theta_t     = PVector -> PVector
minimizer PVector
theta_delta
        zv :: Double
zv          = S -> Array D Ix1 Double -> PVector
forall r e r' ix.
(Manifest r e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
A.computeAs S
A.S ((Double, Array D Ix1 Double) -> Array D Ix1 Double
forall a b. (a, b) -> b
snd ((Double, Array D Ix1 Double) -> Array D Ix1 Double)
-> (Double, Array D Ix1 Double) -> Array D Ix1 Double
forall a b. (a -> b) -> a -> b
$ Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (Double, Array D Ix1 Double)
gradNLL Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta_t) PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix
        zvs :: Array D Ix1 Double
zvs         = (Double, Array D Ix1 Double) -> Array D Ix1 Double
forall a b. (a, b) -> b
snd ((Double, Array D Ix1 Double) -> Array D Ix1 Double)
-> (Double, Array D Ix1 Double) -> Array D Ix1 Double
forall a b. (a -> b) -> a -> b
$ Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (Double, Array D Ix1 Double)
gradNLL Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta_t
        inv_slope' :: Double
inv_slope'  = Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
4.0 (Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0.0625 (Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
forall a. Num a => a -> a
abs (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (Double
tau Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
stdErr_i Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
zv))
        nll_cond :: Double
nll_cond    = Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
nll Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta_t
        acc' :: ([Double], [PVector])
acc'        = if Double
nll_cond Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
nll_opt Bool -> Bool -> Bool
|| ( (Bool -> Bool
not(Bool -> Bool) -> ([Double] -> Bool) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
.[Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null) [Double]
taus Bool -> Bool -> Bool
&& Double
tau Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Double
forall a. HasCallStack => [a] -> a
head [Double]
taus ) Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
tau
                         then ([Double], [PVector])
acc
                         else (Double
tauDouble -> [Double] -> [Double]
forall a. a -> [a] -> [a]
:[Double]
taus, PVector
theta_tPVector -> [PVector] -> [PVector]
forall a. a -> [a] -> [a]
:[PVector]
thetas)
        tau :: Double
tau         = Double -> Double
forall a. Num a => a -> a
signum Double
delta Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
nll_cond Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
nll_opt)

-- Based on https://insysbio.github.io/LikelihoodProfiler.jl/latest/
-- Borisov, Ivan, and Evgeny Metelkin. "Confidence intervals by constrained optimization—An algorithm and software package for practical identifiability analysis in systems biology." PLOS Computational Biology 16.12 (2020): e1008495.
getProfileCnstr :: Distribution
                -> Maybe Double
                -> SRMatrix
                -> PVector
                -> Fix SRTree
                -> PVector
                -> Double -> Double
                -> Int
                -> Either PVector ProfileT
getProfileCnstr :: Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
-> Double
-> Ix1
-> Either PVector ProfileT
getProfileCnstr Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta Double
stdErr_i Double
tau_max Ix1
ix
  | Double
stdErr_i Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0.0 = ProfileT -> Either PVector ProfileT
forall a. a -> Either PVector a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (ProfileT -> Either PVector ProfileT)
-> ProfileT -> Either PVector ProfileT
forall a b. (a -> b) -> a -> b
$ PVector
-> SRMatrix
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> ProfileT
ProfileT PVector
taus SRMatrix
thetas Double
theta_i (Double -> Double -> Double
forall a b. a -> b -> a
const Double
theta_i) (Double -> Double -> Double
forall a b. a -> b -> a
const Double
tau_max)
  | Bool
otherwise       = ProfileT -> Either PVector ProfileT
forall a. a -> Either PVector a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (ProfileT -> Either PVector ProfileT)
-> ProfileT -> Either PVector ProfileT
forall a b. (a -> b) -> a -> b
$ PVector
-> SRMatrix
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> ProfileT
ProfileT PVector
taus SRMatrix
thetas Double
theta_i Double -> Double
forall {a}. (Ord a, Num a) => a -> Double
tau2theta (Double -> Double -> Double
forall a b. a -> b -> a
const Double
tau_max)
  where
    taus :: PVector
taus     = Comp -> [Double] -> PVector
forall r e. Manifest r e => Comp -> [e] -> Vector r e
A.fromList Comp
compMode [-Double
tau_max, Double
tau_max]
    theta' :: [Double]
theta'   = PVector -> [Double]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList PVector
theta
    thetas :: SRMatrix
thetas   = Comp -> [ListItem Ix2 Double] -> SRMatrix
forall r ix e.
(HasCallStack, Ragged L ix e, Manifest r e) =>
Comp -> [ListItem ix e] -> Array r ix e
A.fromLists' Comp
compMode [[Double]
ListItem Ix2 Double
theta', [Double]
ListItem Ix2 Double
theta']
    theta_i :: Double
theta_i  = PVector
theta PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix
    getPoint :: Bool -> Double
getPoint = Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
-> Ix1
-> Bool
-> Double
getEndPoint Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta Double
tau_max Ix1
ix
    leftPt :: Double
leftPt   = Bool -> Double
getPoint Bool
True
    rightPt :: Double
rightPt  = Bool -> Double
getPoint Bool
False
    tau2theta :: a -> Double
tau2theta a
tau = if a
tau a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 then Double
leftPt else Double
rightPt

getEndPoint :: Distribution -> Maybe Double -> A.Array A.S Ix2 Double -> A.Array A.S A.Ix1 Double -> Fix SRTree -> A.Array A.S A.Ix1 Double -> Double -> Int -> Bool -> Double
getEndPoint :: Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
-> Ix1
-> Bool
-> Double
getEndPoint Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta Double
tau_max Ix1
ix Bool
isLeft =
  case AugLagProblem -> Vector Double -> Either Result Solution
minimizeAugLag AugLagProblem
problem (PVector -> Vector Double
forall ix e. Index ix => Array S ix e -> Vector e
A.toStorableVector PVector
theta_opt) of
            Right Solution
sol -> Solution -> Vector Double
solutionParams Solution
sol Vector Double -> Ix1 -> Double
forall a. Storable a => Vector a -> Ix1 -> a
VS.! Ix1
ix
            Left Result
e    -> Result -> Double -> Double
forall a b. Show a => a -> b -> b
traceShow Result
e (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ PVector
theta_opt PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix
  where
    (A.Sz1 Ix1
n) = PVector -> Sz Ix1
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
A.size PVector
theta

    theta_opt :: PVector
theta_opt = (PVector, Double) -> PVector
forall a b. (a, b) -> a
fst ((PVector, Double) -> PVector) -> (PVector, Double) -> PVector
forall a b. (a -> b) -> a -> b
$ Distribution
-> Maybe Double
-> Ix1
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (PVector, Double)
minimizeNLLNonUnique Distribution
dist Maybe Double
mSErr Ix1
100 SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta
    nll_opt :: Double
nll_opt   = Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
nll Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta_opt
    loss_crit :: Double
loss_crit = Double
nll_opt Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
tau_max

    loss :: Vector Double -> Double
loss      = Double -> Double -> Double
forall a. Num a => a -> a -> a
subtract Double
loss_crit (Double -> Double)
-> (Vector Double -> Double) -> Vector Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
nll Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree (PVector -> Double)
-> (Vector Double -> PVector) -> Vector Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Comp -> Vector Double -> PVector
forall e. Comp -> Vector e -> Vector S e
A.fromStorableVector Comp
compMode
    obj :: Vector Double -> Double
obj       = (if Bool
isLeft then Double -> Double
forall a. a -> a
id else Double -> Double
forall a. Num a => a -> a
negate) (Double -> Double)
-> (Vector Double -> Double) -> Vector Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector Double -> Ix1 -> Double
forall a. Storable a => Vector a -> Ix1 -> a
VS.! Ix1
ix)

    stop :: NonEmpty StoppingCondition
stop       = Double -> StoppingCondition
ObjectiveRelativeTolerance Double
1e-4 StoppingCondition
-> [StoppingCondition] -> NonEmpty StoppingCondition
forall a. a -> [a] -> NonEmpty a
:| []
    localAlg :: LocalAlgorithm
localAlg   = (Vector Double -> Double)
-> [Bounds] -> Maybe InitialStep -> LocalAlgorithm
NELDERMEAD Vector Double -> Double
obj [] Maybe InitialStep
forall a. Maybe a
Nothing
    local :: LocalProblem
local      = Word
-> NonEmpty StoppingCondition -> LocalAlgorithm -> LocalProblem
LocalProblem (Ix1 -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral Ix1
n) NonEmpty StoppingCondition
stop LocalAlgorithm
localAlg
    constraint :: InequalityConstraint (Vector Double -> Double) v
constraint = Constraint (Vector Double -> Double) v
-> Double -> InequalityConstraint (Vector Double -> Double) v
forall s v. Constraint s v -> Double -> InequalityConstraint s v
InequalityConstraint ((Vector Double -> Double) -> Constraint (Vector Double -> Double) v
forall s v. s -> Constraint s v
Scalar Vector Double -> Double
loss) Double
1e-6

    problem :: AugLagProblem
problem = EqualityConstraints
-> EqualityConstraintsD -> AugLagAlgorithm -> AugLagProblem
AugLagProblem [] [] (LocalProblem
-> InequalityConstraints
-> InequalityConstraintsD
-> AugLagAlgorithm
AUGLAG_LOCAL LocalProblem
local [InequalityConstraint (Vector Double -> Double) VectorConstraint
forall {v}. InequalityConstraint (Vector Double -> Double) v
constraint] [])
{-# INLINE getEndPoint #-}

-- Based on
-- Jian-Shen Chen & Robert I Jennrich (2002) Simple Accurate Approximation of Likelihood Profiles,
-- Journal of Computational and Graphical Statistics, 11:3, 714-732, DOI: 10.1198/106186002493
getProfileODE :: Distribution
           -> Maybe Double
           -> SRMatrix
           -> PVector
           -> Fix SRTree
           -> PVector
           -> Double
           -> CI
           -> Double
           -> Int
           -> Either PVector ProfileT
getProfileODE :: Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
-> CI
-> Double
-> Ix1
-> Either PVector ProfileT
getProfileODE Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta Double
stdErr_i CI
estCI Double
tau_max Ix1
ix
  | Double
stdErr_i Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0.0 = ProfileT -> Either PVector ProfileT
forall a. a -> Either PVector a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ProfileT
dflt
  | Bool
otherwise = let (Comp -> [Double] -> PVector
forall r e. Manifest r e => Comp -> [e] -> Vector r e
A.fromList Comp
compMode -> PVector
taus, Comp -> [ListItem Ix2 Double] -> SRMatrix
forall r ix e.
(HasCallStack, Ragged L ix e, Manifest r e) =>
Comp -> [ListItem ix e] -> Array r ix e
A.fromLists' Comp
compMode ([[Double]] -> SRMatrix)
-> ([PVector] -> [[Double]]) -> [PVector] -> SRMatrix
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (PVector -> [Double]) -> [PVector] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map PVector -> [Double]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList -> SRMatrix
thetas) = ([Double], [PVector])
solLeft ([Double], [PVector])
-> ([Double], [PVector]) -> ([Double], [PVector])
forall a. Semigroup a => a -> a -> a
<> ([Double
0], [PVector
theta_opt]) ([Double], [PVector])
-> ([Double], [PVector]) -> ([Double], [PVector])
forall a. Semigroup a => a -> a -> a
<> ([Double], [PVector])
solRight
                    (Double -> Double
tau2theta, Double -> Double
theta2tau) = PVector
-> SRMatrix
-> Double
-> Double
-> Ix1
-> (Double -> Double, Double -> Double)
createSplines PVector
taus SRMatrix
thetas Double
stdErr_i Double
tau_max Ix1
ix
                in ProfileT -> Either PVector ProfileT
forall a. a -> Either PVector a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (ProfileT -> Either PVector ProfileT)
-> ProfileT -> Either PVector ProfileT
forall a b. (a -> b) -> a -> b
$ PVector
-> SRMatrix
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> ProfileT
ProfileT PVector
taus SRMatrix
thetas Double
optTh Double -> Double
tau2theta Double -> Double
theta2tau
  where
    dflt :: ProfileT
dflt      = PVector
-> SRMatrix
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> ProfileT
ProfileT (Comp -> [Double] -> PVector
forall r e. Manifest r e => Comp -> [e] -> Vector r e
A.fromList Comp
compMode [-Double
tau_max, Double
tau_max]) (Comp -> [ListItem Ix2 Double] -> SRMatrix
forall r ix e.
(HasCallStack, Ragged L ix e, Manifest r e) =>
Comp -> [ListItem ix e] -> Array r ix e
A.fromLists' Comp
compMode [[Double]
ListItem Ix2 Double
theta', [Double]
ListItem Ix2 Double
theta']) (PVector
theta PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix) (Double -> Double -> Double
forall a b. a -> b -> a
const (PVector
theta PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix)) (Double -> Double -> Double
forall a b. a -> b -> a
const Double
tau_max)
    minimizer :: PVector -> PVector
minimizer = (PVector, Double) -> PVector
forall a b. (a, b) -> a
fst ((PVector, Double) -> PVector)
-> (PVector -> (PVector, Double)) -> PVector -> PVector
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Distribution
-> Maybe Double
-> Ix1
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (PVector, Double)
minimizeNLLNonUnique Distribution
dist Maybe Double
mSErr Ix1
100 SRMatrix
xss PVector
ys Fix SRTree
tree
    grader :: PVector -> Array D Ix1 Double
grader    = (Double, Array D Ix1 Double) -> Array D Ix1 Double
forall a b. (a, b) -> b
snd ((Double, Array D Ix1 Double) -> Array D Ix1 Double)
-> (PVector -> (Double, Array D Ix1 Double))
-> PVector
-> Array D Ix1 Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (Double, Array D Ix1 Double)
gradNLLNonUnique Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree
    theta_opt :: PVector
theta_opt = PVector -> PVector
minimizer PVector
theta
    theta' :: [Double]
theta'    = PVector -> [Double]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList PVector
theta
    nll_opt :: Double
nll_opt   = Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
nll Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta_opt
    optTh :: Double
optTh     = PVector
theta_opt PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix
    p' :: Ix1
p'        = Ix1
pIx1 -> Ix1 -> Ix1
forall a. Num a => a -> a -> a
+Ix1
1
    (A.Sz1 Ix1
p) = PVector -> Sz Ix1
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
A.size PVector
theta
    sErr :: Double
sErr      = Double -> Maybe Double -> Double
forall a. a -> Maybe a -> a
fromMaybe Double
1 Maybe Double
mSErr
    getHess :: PVector -> SRMatrix
getHess   = Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> SRMatrix
hessianNLL Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree

    odeFun :: Double -> p -> PVector -> PVector
odeFun Double
gamma p
_ PVector
u =
        let grad :: Array D Ix1 Double
grad     = PVector -> Array D Ix1 Double
grader PVector
u
            w :: SRMatrix
w        = Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> SRMatrix
hessianNLL Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
u
            m :: SRMatrix
m        = Comp -> Sz Ix2 -> (Ix2 -> Double) -> SRMatrix
forall r ix e.
Load r ix e =>
Comp -> Sz ix -> (ix -> e) -> Array r ix e
A.makeArray Comp
compMode (Ix2 -> Sz Ix2
forall ix. Index ix => ix -> Sz ix
A.Sz (Ix1
p' Ix1 -> Ix1 -> Ix2
:. Ix1
p'))
                         (\ (Ix1
i :. Ix1
j) -> if | Ix1
iIx1 -> Ix1 -> Bool
forall a. Ord a => a -> a -> Bool
<Ix1
p Bool -> Bool -> Bool
&& Ix1
jIx1 -> Ix1 -> Bool
forall a. Ord a => a -> a -> Bool
<Ix1
p -> SRMatrix
w SRMatrix -> Ix2 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! (Ix1
i Ix1 -> Ix1 -> Ix2
:. Ix1
j)
                                           | Ix1
iIx1 -> Ix1 -> Bool
forall a. Eq a => a -> a -> Bool
==Ix1
ix      -> Double
1
                                           | Ix1
jIx1 -> Ix1 -> Bool
forall a. Eq a => a -> a -> Bool
==Ix1
ix      -> Double
1
                                           | Bool
otherwise  -> Double
0
                         )

            v :: PVector
v        = S -> Array DL Ix1 Double -> PVector
forall r e r' ix.
(Manifest r e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
A.computeAs S
A.S (Array DL Ix1 Double -> PVector) -> Array DL Ix1 Double -> PVector
forall a b. (a -> b) -> a -> b
$ Array D Ix1 Double -> Double -> Array DL Ix1 Double
forall r e.
(Size r, Load r Ix1 e) =>
Vector r e -> e -> Vector DL e
A.snoc ((Double -> Double) -> Array D Ix1 Double -> Array D Ix1 Double
forall ix r e' e.
(Index ix, Source r e') =>
(e' -> e) -> Array r ix e' -> Array D ix e
A.map (Double -> Double -> Double
forall a. Num a => a -> a -> a
*(-Double
gamma)) Array D Ix1 Double
grad) Double
1
            dotTheta :: PVector
dotTheta = IO PVector -> PVector
forall a. IO a -> a
unsafePerformIO (IO PVector -> PVector) -> IO PVector -> PVector
forall a b. (a -> b) -> a -> b
$ SRMatrix -> PVector -> IO PVector
forall (m :: * -> *).
(PrimMonad m, MonadThrow m, MonadIO m) =>
SRMatrix -> PVector -> m PVector
luSolve SRMatrix
m PVector
v
        in Comp -> Vector Double -> PVector
forall e. Comp -> Vector e -> Vector S e
A.fromStorableVector Comp
compMode (Vector Double -> PVector) -> Vector Double -> PVector
forall a b. (a -> b) -> a -> b
$ Vector Double -> Vector Double
forall a. Storable a => Vector a -> Vector a
VS.init (Vector Double -> Vector Double) -> Vector Double -> Vector Double
forall a b. (a -> b) -> a -> b
$ PVector -> Vector Double
forall ix e. Index ix => Array S ix e -> Vector e
A.toStorableVector PVector
dotTheta
    tsHi :: [Double]
tsHi = Ix1 -> (Double, Double) -> [Double]
linSpace Ix1
50 (Double
optTh, CI -> Double
upper_ CI
estCI)
    tsLo :: [Double]
tsLo = Ix1 -> (Double, Double) -> [Double]
linSpace Ix1
50 (Double
optTh, CI -> Double
lower_ CI
estCI)
    scanOn :: Double -> [Double] -> ([Double], [PVector])
scanOn Double
sig = ((Double, PVector) -> ([Double], [PVector]))
-> [(Double, PVector)] -> ([Double], [PVector])
forall m a. Monoid m => (a -> m) -> [a] -> m
forall (t :: * -> *) m a.
(Foldable t, Monoid m) =>
(a -> m) -> t a -> m
foldMap (Double -> (Double, PVector) -> ([Double], [PVector])
forall {p} {a}. p -> (a, PVector) -> ([Double], [PVector])
calcTau Double
sig) ([(Double, PVector)] -> ([Double], [PVector]))
-> ([Double] -> [(Double, PVector)])
-> [Double]
-> ([Double], [PVector])
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(Double, PVector)] -> [(Double, PVector)]
forall a. [a] -> [a]
f ([(Double, PVector)] -> [(Double, PVector)])
-> ([Double] -> [(Double, PVector)])
-> [Double]
-> [(Double, PVector)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Double, PVector) -> Double -> (Double, PVector))
-> (Double, PVector) -> [Double] -> [(Double, PVector)]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl ((Double -> PVector -> PVector)
-> (Double, PVector) -> Double -> (Double, PVector)
rk (Double -> Double -> PVector -> PVector
forall {p}. Double -> p -> PVector -> PVector
odeFun Double
sig)) (Double
optTh, PVector
theta_opt)
                    where f :: [a] -> [a]
f = if Double
sigDouble -> Double -> Bool
forall a. Eq a => a -> a -> Bool
==Double
1 then [a] -> [a]
forall a. a -> a
id else [a] -> [a]
forall a. [a] -> [a]
reverse
    solRight :: ([Double], [PVector])
solRight = Double -> [Double] -> ([Double], [PVector])
scanOn Double
1 [Double]
tsHi
    solLeft :: ([Double], [PVector])
solLeft  = Double -> [Double] -> ([Double], [PVector])
scanOn (-Double
1) [Double]
tsLo
    calcTau :: p -> (a, PVector) -> ([Double], [PVector])
calcTau p
s (a, PVector)
t = let nll_i :: Double
nll_i = Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> Double
nll Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree (PVector -> Double) -> PVector -> Double
forall a b. (a -> b) -> a -> b
$ (a, PVector) -> PVector
forall a b. (a, b) -> b
snd (a, PVector)
t
                      z :: Double
z     = Double -> Double
forall a. Num a => a -> a
signum (((a, PVector) -> PVector
forall a b. (a, b) -> b
snd (a, PVector)
t PVector -> Ix1 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
optTh) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
nll_i Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
nll_opt)
                   in if Double
z Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
z then ([], []) else ([Double
z], [(a, PVector) -> PVector
forall a b. (a, b) -> b
snd (a, PVector)
t])

rk :: (Double -> PVector -> PVector) -> (Double, PVector) -> Double -> (Double, PVector)
rk :: (Double -> PVector -> PVector)
-> (Double, PVector) -> Double -> (Double, PVector)
rk Double -> PVector -> PVector
f (Double
t, PVector
y) Double
t' = (Double
t', PVector
y PVector -> PVector -> PVector
forall ix r e.
(HasCallStack, Index ix, Numeric r e) =>
Array r ix e -> Array r ix e -> Array r ix e
!+! ((Double
1.0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
6.0) Double -> PVector -> PVector
forall ix r e.
(Index ix, Numeric r e) =>
e -> Array r ix e -> Array r ix e
*. PVector
h' PVector -> PVector -> PVector
forall ix r e.
(Index ix, Numeric r e) =>
Array r ix e -> Array r ix e -> Array r ix e
!*! (PVector
k1 PVector -> PVector -> PVector
forall ix r e.
(HasCallStack, Index ix, Numeric r e) =>
Array r ix e -> Array r ix e -> Array r ix e
!+! (Double
2.0 Double -> PVector -> PVector
forall ix r e.
(Index ix, Numeric r e) =>
e -> Array r ix e -> Array r ix e
*. PVector
k2) PVector -> PVector -> PVector
forall ix r e.
(HasCallStack, Index ix, Numeric r e) =>
Array r ix e -> Array r ix e -> Array r ix e
!+! (Double
2.0 Double -> PVector -> PVector
forall ix r e.
(Index ix, Numeric r e) =>
e -> Array r ix e -> Array r ix e
*. PVector
k3) PVector -> PVector -> PVector
forall ix r e.
(HasCallStack, Index ix, Numeric r e) =>
Array r ix e -> Array r ix e -> Array r ix e
!+! PVector
k4)))
  where
    h :: Double
h  = Double
t' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t
    h', k1, k2, k3, k4 :: PVector
    h' :: PVector
h' = Comp -> Sz Ix1 -> Double -> PVector
forall r ix e. Load r ix e => Comp -> Sz ix -> e -> Array r ix e
A.replicate Comp
compMode (PVector -> Sz Ix1
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
A.size PVector
y) Double
h
    k1 :: PVector
k1 = Double -> PVector -> PVector
f Double
t PVector
y
    k2 :: PVector
k2 = Double -> PVector -> PVector
f (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.5Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
h) (S -> Array D Ix1 Double -> PVector
forall r e r' ix.
(Manifest r e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
A.computeAs S
A.S (Array D Ix1 Double -> PVector) -> Array D Ix1 Double -> PVector
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double -> Double)
-> PVector -> PVector -> PVector -> Array D Ix1 Double
forall ix r1 e1 r2 e2 r3 e3 e.
(Index ix, Source r1 e1, Source r2 e2, Source r3 e3) =>
(e1 -> e2 -> e3 -> e)
-> Array r1 ix e1
-> Array r2 ix e2
-> Array r3 ix e3
-> Array D ix e
A.zipWith3 (Double -> Double -> Double -> Double -> Double
forall {a}. Num a => a -> a -> a -> a -> a
g Double
0.5) PVector
y PVector
h' PVector
k1) -- (y !+! 0.5*.h' A.!*! k1)
    k3 :: PVector
k3 = Double -> PVector -> PVector
f (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.5Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
h) (S -> Array D Ix1 Double -> PVector
forall r e r' ix.
(Manifest r e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
A.computeAs S
A.S (Array D Ix1 Double -> PVector) -> Array D Ix1 Double -> PVector
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double -> Double)
-> PVector -> PVector -> PVector -> Array D Ix1 Double
forall ix r1 e1 r2 e2 r3 e3 e.
(Index ix, Source r1 e1, Source r2 e2, Source r3 e3) =>
(e1 -> e2 -> e3 -> e)
-> Array r1 ix e1
-> Array r2 ix e2
-> Array r3 ix e3
-> Array D ix e
A.zipWith3 (Double -> Double -> Double -> Double -> Double
forall {a}. Num a => a -> a -> a -> a -> a
g Double
0.5) PVector
y PVector
h' PVector
k2) -- (y !+! 0.5*.h' A.!*! k2)
    k4 :: PVector
k4 = Double -> PVector -> PVector
f (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1.0Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
h) (S -> Array D Ix1 Double -> PVector
forall r e r' ix.
(Manifest r e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
A.computeAs S
A.S (Array D Ix1 Double -> PVector) -> Array D Ix1 Double -> PVector
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double -> Double)
-> PVector -> PVector -> PVector -> Array D Ix1 Double
forall ix r1 e1 r2 e2 r3 e3 e.
(Index ix, Source r1 e1, Source r2 e2, Source r3 e3) =>
(e1 -> e2 -> e3 -> e)
-> Array r1 ix e1
-> Array r2 ix e2
-> Array r3 ix e3
-> Array D ix e
A.zipWith3 (Double -> Double -> Double -> Double -> Double
forall {a}. Num a => a -> a -> a -> a -> a
g Double
1.0) PVector
y PVector
h' PVector
k3) -- (y !+! 1.0*.h'!*!k3)
    g :: a -> a -> a -> a -> a
g a
a a
yi a
hi a
ki = a
yi a -> a -> a
forall a. Num a => a -> a -> a
+ a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
hi a -> a -> a
forall a. Num a => a -> a -> a
* a
ki
{-# INLINE rk #-}

-- tau0, tau1  theta0, thetaX = tau1 theta0 / tau0
getStatsFromModel :: Distribution -> Maybe Double -> SRMatrix -> PVector -> Fix SRTree -> PVector -> BasicStats
getStatsFromModel :: Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> BasicStats
getStatsFromModel Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta = SRMatrix -> SRMatrix -> PVector -> BasicStats
MkStats SRMatrix
cov SRMatrix
corr PVector
stdErr
  where
    (A.Sz1 Ix1
k) = PVector -> Sz Ix1
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
A.size PVector
theta
    (A.Sz1 Ix1
n) = PVector -> Sz Ix1
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
A.size PVector
ys
    nParams :: Sz Ix1
nParams = Ix1 -> Sz Ix1
forall a b. (Integral a, Num b) => a -> b
fromIntegral Ix1
k
    ssr :: Double
ssr  = SRMatrix -> PVector -> Fix SRTree -> PVector -> Double
sse SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta
    ident :: SRMatrix
ident = S -> Array DL Ix2 Double -> SRMatrix
forall r e r' ix.
(Manifest r e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
A.computeAs S
A.S (Array DL Ix2 Double -> SRMatrix)
-> Array DL Ix2 Double -> SRMatrix
forall a b. (a -> b) -> a -> b
$ Sz Ix1 -> Array DL Ix2 Double
forall e. Num e => Sz Ix1 -> Matrix DL e
identityMatrix Sz Ix1
nParams

    -- only for gaussian
    sErr :: Double
sErr  = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
ssr Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Ix1 -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Ix1
n Ix1 -> Ix1 -> Ix1
forall a. Num a => a -> a -> a
- Ix1
k)

    hess :: SRMatrix
hess    = Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> SRMatrix
hessianNLL Distribution
dist Maybe Double
mSErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
theta
    -- cov     = catch (unsafePerformIO (invChol hess)) (\e -> trace "cov NegDef" $ pure ident)
    fexcept :: (A.PrimMonad m, A.MonadThrow m, A.MonadIO m) => A.SomeException -> m SRMatrix
    fexcept :: forall (m :: * -> *).
(PrimMonad m, MonadThrow m, MonadIO m) =>
SomeException -> m SRMatrix
fexcept SomeException
e = String -> m SRMatrix -> m SRMatrix
forall a. String -> a -> a
trace String
"cov NegDef" (m SRMatrix -> m SRMatrix) -> m SRMatrix -> m SRMatrix
forall a b. (a -> b) -> a -> b
$ SRMatrix -> m SRMatrix
forall a. a -> m a
forall (f :: * -> *) a. Applicative f => a -> f a
pure SRMatrix
ident
    cov :: SRMatrix
cov     = IO SRMatrix -> SRMatrix
forall a. IO a -> a
unsafePerformIO (IO SRMatrix -> SRMatrix) -> IO SRMatrix -> SRMatrix
forall a b. (a -> b) -> a -> b
$ IO SRMatrix -> (SomeException -> IO SRMatrix) -> IO SRMatrix
forall e a.
(HasCallStack, Exception e) =>
IO a -> (e -> IO a) -> IO a
forall (m :: * -> *) e a.
(MonadCatch m, HasCallStack, Exception e) =>
m a -> (e -> m a) -> m a
catch (SRMatrix -> IO SRMatrix
forall (m :: * -> *).
(PrimMonad m, MonadThrow m, MonadIO m) =>
SRMatrix -> m SRMatrix
invChol SRMatrix
hess) SomeException -> IO SRMatrix
forall (m :: * -> *).
(PrimMonad m, MonadThrow m, MonadIO m) =>
SomeException -> m SRMatrix
fexcept

    stdErr :: PVector
stdErr   = Comp -> Sz Ix1 -> (Ix1 -> Double) -> PVector
forall r ix e.
Load r ix e =>
Comp -> Sz ix -> (ix -> e) -> Array r ix e
A.makeArray Comp
compMode (Ix1 -> Sz Ix1
A.Sz1 Ix1
k) (\Ix1
ix -> Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ SRMatrix
cov SRMatrix -> Ix2 -> Double
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! (Ix1
ix Ix1 -> Ix1 -> Ix2
:. Ix1
ix))
    stdErrSq :: SRMatrix
stdErrSq = case PVector -> PVector -> Either SomeException SRMatrix
forall (m :: * -> *).
MonadThrow m =>
PVector -> PVector -> m SRMatrix
outer PVector
stdErr PVector
stdErr of
                 Left SomeException
_  -> String -> SRMatrix
forall a. HasCallStack => String -> a
error String
"stdErr size mismatch?"
                 Right SRMatrix
v -> SRMatrix
v

    corr :: SRMatrix
corr     = S -> Array D Ix2 Double -> SRMatrix
forall r e r' ix.
(Manifest r e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
A.computeAs S
A.S (Array D Ix2 Double -> SRMatrix) -> Array D Ix2 Double -> SRMatrix
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double)
-> SRMatrix -> SRMatrix -> Array D Ix2 Double
forall ix r1 e1 r2 e2 e.
(Index ix, Source r1 e1, Source r2 e2) =>
(e1 -> e2 -> e) -> Array r1 ix e1 -> Array r2 ix e2 -> Array D ix e
A.zipWith Double -> Double -> Double
forall a. Fractional a => a -> a -> a
(/) SRMatrix
cov SRMatrix
stdErrSq

-- Create splines for profile-t
createSplines :: PVector -> SRMatrix -> Double -> Double -> Int -> (Double -> Double, Double -> Double)
createSplines :: PVector
-> SRMatrix
-> Double
-> Double
-> Ix1
-> (Double -> Double, Double -> Double)
createSplines PVector
taus SRMatrix
thetas Double
se Double
tau_max Ix1
ix
  | Ix1
n Ix1 -> Ix1 -> Bool
forall a. Ord a => a -> a -> Bool
< Ix1
2     = ([(Double, Double)] -> Double -> Double
genSplineFun [(-Double
tau_max, -Double
se), (Double
tau_max, Double
se)], [(Double, Double)] -> Double -> Double
genSplineFun [(-Double
se, Double
0), (Double
se, Double
1)])
  | Bool
otherwise = (Double -> Double
tau2theta, Double -> Double
theta2tau)
  where
    (A.Sz Ix1
n)   = PVector -> Sz Ix1
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
A.size PVector
taus
    cols :: PVector
cols       = Ix1 -> SRMatrix -> PVector
getCol Ix1
ix SRMatrix
thetas
    nubOnFirst :: [(Double, b)] -> [(Double, b)]
nubOnFirst = ((Double, b) -> (Double, b) -> Bool)
-> [(Double, b)] -> [(Double, b)]
forall a. (a -> a -> Bool) -> [a] -> [a]
nubBy (\(Double, b)
x (Double, b)
y -> (Double, b) -> Double
forall a b. (a, b) -> a
fst (Double, b)
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== (Double, b) -> Double
forall a b. (a, b) -> a
fst (Double, b)
y)
    tau2theta :: Double -> Double
tau2theta  = [(Double, Double)] -> Double -> Double
genSplineFun ([(Double, Double)] -> Double -> Double)
-> [(Double, Double)] -> Double -> Double
forall a b. (a -> b) -> a -> b
$ [(Double, Double)] -> [(Double, Double)]
forall {b}. [(Double, b)] -> [(Double, b)]
nubOnFirst ([(Double, Double)] -> [(Double, Double)])
-> [(Double, Double)] -> [(Double, Double)]
forall a b. (a -> b) -> a -> b
$ PVector -> PVector -> [(Double, Double)]
sortOnFirst PVector
taus PVector
cols
    theta2tau :: Double -> Double
theta2tau  = [(Double, Double)] -> Double -> Double
genSplineFun ([(Double, Double)] -> Double -> Double)
-> [(Double, Double)] -> Double -> Double
forall a b. (a -> b) -> a -> b
$ [(Double, Double)] -> [(Double, Double)]
forall {b}. [(Double, b)] -> [(Double, b)]
nubOnFirst ([(Double, Double)] -> [(Double, Double)])
-> [(Double, Double)] -> [(Double, Double)]
forall a b. (a -> b) -> a -> b
$ PVector -> PVector -> [(Double, Double)]
sortOnFirst PVector
cols PVector
taus

getCol :: Int -> SRMatrix -> PVector
getCol :: Ix1 -> SRMatrix -> PVector
getCol Ix1
ix SRMatrix
mtx = SRMatrix -> Array B Ix1 PVector
getCols SRMatrix
mtx Array B Ix1 PVector -> Ix1 -> PVector
forall r ix e.
(HasCallStack, Manifest r e, Index ix) =>
Array r ix e -> ix -> e
A.! Ix1
ix
{-# inline getCol #-}

sortOnFirst :: PVector -> PVector -> [(Double, Double)]
sortOnFirst :: PVector -> PVector -> [(Double, Double)]
sortOnFirst PVector
xs PVector
ys = ((Double, Double) -> Double)
-> [(Double, Double)] -> [(Double, Double)]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (Double, Double) -> Double
forall a b. (a, b) -> a
fst ([(Double, Double)] -> [(Double, Double)])
-> [(Double, Double)] -> [(Double, Double)]
forall a b. (a -> b) -> a -> b
$ [Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip (PVector -> [Double]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList PVector
xs) (PVector -> [Double]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList PVector
ys)
{-# inline sortOnFirst #-}

splinesSketches :: Double -> PVector -> PVector -> (Double -> Double) -> (Double -> Double)
splinesSketches :: Double
-> PVector -> PVector -> (Double -> Double) -> Double -> Double
splinesSketches Double
tauScale (PVector -> [Double]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList -> [Double]
tau) (PVector -> [Double]
forall ix r e. (Index ix, Source r e) => Array r ix e -> [e]
A.toList -> [Double]
theta) Double -> Double
theta2tau
  | [Double] -> Ix1
forall a. [a] -> Ix1
forall (t :: * -> *) a. Foldable t => t a -> Ix1
length [Double]
tau Ix1 -> Ix1 -> Bool
forall a. Ord a => a -> a -> Bool
< Ix1
2 = Double -> Double
forall a. a -> a
id
  | Bool
otherwise      = [(Double, Double)] -> Double -> Double
genSplineFun [(Double, Double)]
gpq
  where
    gpq :: [(Double, Double)]
gpq = ((Double, Double) -> Double)
-> [(Double, Double)] -> [(Double, Double)]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (Double, Double) -> Double
forall a b. (a, b) -> a
fst [(Double
x, Double -> Double
forall a. Floating a => a -> a
acos Double
y') | (Double
x, Double
y) <- [Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
tau [Double]
theta
                                   , let y' :: Double
y' = Double -> Double
theta2tau Double
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tauScale
                                   , Double -> Double
forall a. Num a => a -> a
abs Double
y' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1 ]

approximateContour :: Int -> Int -> [ProfileT] -> Int -> Int -> Double -> [(Double, Double)]
approximateContour :: Ix1
-> Ix1 -> [ProfileT] -> Ix1 -> Ix1 -> Double -> [(Double, Double)]
approximateContour Ix1
nParams Ix1
nPoints [ProfileT]
profs Ix1
ix1 Ix1
ix2 Double
alpha = Double -> [(Double, Double)]
go Double
0
  where
    -- get the info for ix1 and ix2
    (ProfileT
prof1, ProfileT
prof2)           = ([ProfileT]
profs [ProfileT] -> Ix1 -> ProfileT
forall a. HasCallStack => [a] -> Ix1 -> a
!! Ix1
ix1, [ProfileT]
profs [ProfileT] -> Ix1 -> ProfileT
forall a. HasCallStack => [a] -> Ix1 -> a
!! Ix1
ix2)
    (Double -> Double
tau2theta1, Double -> Double
theta2tau1) = (ProfileT -> Double -> Double
_tau2theta ProfileT
prof1, ProfileT -> Double -> Double
_theta2tau ProfileT
prof1)
    (Double -> Double
tau2theta2, Double -> Double
theta2tau2) = (ProfileT -> Double -> Double
_tau2theta ProfileT
prof2, ProfileT -> Double -> Double
_theta2tau ProfileT
prof2)

    -- calculate the spline for A-D
    tauScale :: Double
tauScale = Double -> Double
forall a. Floating a => a -> a
sqrt (Ix1 -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Ix1
nParams Double -> Double -> Double
forall a. Num a => a -> a -> a
* FDistribution -> Double -> Double
forall d. ContDistr d => d -> Double -> Double
quantile (Ix1 -> Ix1 -> FDistribution
fDistribution Ix1
nParams (Ix1
nPoints Ix1 -> Ix1 -> Ix1
forall a. Num a => a -> a -> a
- Ix1
nParams)) (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha))
    splineG1 :: Double -> Double
splineG1 = Double
-> PVector -> PVector -> (Double -> Double) -> Double -> Double
splinesSketches Double
tauScale (ProfileT -> PVector
_taus ProfileT
prof1) (Ix1 -> SRMatrix -> PVector
getCol Ix1
ix2 (ProfileT -> SRMatrix
_thetas ProfileT
prof1)) Double -> Double
theta2tau2
    splineG2 :: Double -> Double
splineG2 = Double
-> PVector -> PVector -> (Double -> Double) -> Double -> Double
splinesSketches Double
tauScale (ProfileT -> PVector
_taus ProfileT
prof2) (Ix1 -> SRMatrix -> PVector
getCol Ix1
ix1 (ProfileT -> SRMatrix
_thetas ProfileT
prof2)) Double -> Double
theta2tau1
    angles :: [(Double, Double)]
angles   = [ (Double
0, Double -> Double
splineG1 Double
1), (Double -> Double
splineG2 Double
1, Double
0), (Double
forall a. Floating a => a
pi, Double -> Double
splineG1 (-Double
1)), (Double -> Double
splineG2 (-Double
1), Double
forall a. Floating a => a
pi) ]
    splineAD :: Double -> Double
splineAD = [(Double, Double)] -> Double -> Double
genSplineFun [(Double, Double)]
points

    applyIfNeg :: (a, b) -> (a, b)
applyIfNeg (a
x, b
y) = if b
y b -> b -> Bool
forall a. Ord a => a -> a -> Bool
< b
0 then (-a
x, -b
y) else (a
x ,b
y)
    points :: [(Double, Double)]
points   = ((Double, Double) -> Double)
-> [(Double, Double)] -> [(Double, Double)]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (Double, Double) -> Double
forall a b. (a, b) -> a
fst
             ([(Double, Double)] -> [(Double, Double)])
-> [(Double, Double)] -> [(Double, Double)]
forall a b. (a -> b) -> a -> b
$ [(Double, Double) -> (Double, Double)
forall {b} {a}. (Ord b, Num b, Num a) => (a, b) -> (a, b)
applyIfNeg ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
y)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2, Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y) | (Double
x, Double
y) <- [(Double, Double)]
angles]
            [(Double, Double)] -> [(Double, Double)] -> [(Double, Double)]
forall a. Semigroup a => a -> a -> a
<> (\(Double
x,Double
y) -> [(Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi, Double
y)]) ([(Double, Double)] -> (Double, Double)
forall a. HasCallStack => [a] -> a
head [(Double, Double)]
points)

    -- generate the points of the curve
    go :: Double -> [(Double, Double)]
go Double
100 = []
    go Double
ix  = (Double
p, Double
q) (Double, Double) -> [(Double, Double)] -> [(Double, Double)]
forall a. a -> [a] -> [a]
: Double -> [(Double, Double)]
go (Double
ixDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1)
      where
        ai :: Double
ai = Double
ix Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
99 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
forall a. Floating a => a
pi
        di :: Double
di = Double -> Double
splineAD Double
ai
        taup :: Double
taup = Double -> Double
forall a. Floating a => a -> a
cos (Double
ai Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
di Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tauScale
        tauq :: Double
tauq = Double -> Double
forall a. Floating a => a -> a
cos (Double
ai Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
di Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tauScale
        p :: Double
p = Double -> Double
tau2theta1 Double
taup
        q :: Double
q = Double -> Double
tau2theta2 Double
tauq