{-# LANGUAGE BangPatterns #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Algorithm.SRTree.Opt 
-- 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.Opt
    where

import Algorithm.SRTree.Likelihoods
import Algorithm.SRTree.NonlinearOpt
import Data.Bifunctor (bimap, second)
import Data.Massiv.Array
import Data.SRTree (Fix (..), SRTree (..), floatConstsToParam, relabelParams, countNodes)
import Data.SRTree.Eval (evalTree, compMode)
import qualified Data.Vector.Storable as VS
import qualified Data.IntMap.Strict as IntMap
import Data.SRTree.Recursion

import Debug.Trace

tree2arr :: Fix SRTree -> IntMap.IntMap (Int, Int, Int, Double)
tree2arr :: Fix SRTree -> IntMap (Int, Int, Int, Double)
tree2arr Fix SRTree
tree = [(Int, (Int, Int, Int, Double))] -> IntMap (Int, Int, Int, Double)
forall a. [(Int, a)] -> IntMap a
IntMap.fromList [(Int, (Int, Int, Int, Double))]
listTree
  where
    height :: Fix SRTree -> Integer
height = (SRTree Integer -> Integer) -> Fix SRTree -> Integer
forall (f :: * -> *) a. Functor f => (f a -> a) -> Fix f -> a
cata SRTree Integer -> Integer
forall {a}. (Num a, Ord a) => SRTree a -> a
alg
      where
        alg :: SRTree a -> a
alg (Var Int
ix) = a
1
        alg (Const Double
x) = a
1
        alg (Param Int
ix) = a
1
        alg (Uni Function
_ a
t) = a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
t
        alg (Bin Op
_ a
l a
r) = a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a -> a
forall a. Ord a => a -> a -> a
max a
l a
r
    listTree :: [(Int, (Int, Int, Int, Double))]
listTree = (forall x. SRTree x -> Int -> SRTree (x, Int))
-> (SRTree [(Int, (Int, Int, Int, Double))]
    -> Int -> [(Int, (Int, Int, Int, Double))])
-> Fix SRTree
-> Int
-> [(Int, (Int, Int, Int, Double))]
forall (f :: * -> *) p a.
Functor f =>
(forall x. f x -> p -> f (x, p))
-> (f a -> p -> a) -> Fix f -> p -> a
accu SRTree x -> Int -> SRTree (x, Int)
forall x. SRTree x -> Int -> SRTree (x, Int)
forall {b} {a}. Num b => SRTree a -> b -> SRTree (a, b)
indexer SRTree [(Int, (Int, Int, Int, Double))]
-> Int -> [(Int, (Int, Int, Int, Double))]
forall {a} {a}.
Num a =>
SRTree [(a, (a, Int, Int, Double))]
-> a -> [(a, (a, Int, Int, Double))]
convert Fix SRTree
tree Int
0

    indexer :: SRTree a -> b -> SRTree (a, b)
indexer (Var Int
ix) b
iy   = Int -> SRTree (a, b)
forall val. Int -> SRTree val
Var Int
ix
    indexer (Const Double
x) b
iy  = Double -> SRTree (a, b)
forall val. Double -> SRTree val
Const Double
x
    indexer (Param Int
ix) b
iy = Int -> SRTree (a, b)
forall val. Int -> SRTree val
Param Int
ix
    indexer (Bin Op
op a
l a
r) b
iy = Op -> (a, b) -> (a, b) -> SRTree (a, b)
forall val. Op -> val -> val -> SRTree val
Bin Op
op (a
l, b
2b -> b -> b
forall a. Num a => a -> a -> a
*b
iyb -> b -> b
forall a. Num a => a -> a -> a
+b
1) (a
r, b
2b -> b -> b
forall a. Num a => a -> a -> a
*b
iyb -> b -> b
forall a. Num a => a -> a -> a
+b
2)
    indexer (Uni Function
f a
t) b
iy = Function -> (a, b) -> SRTree (a, b)
forall val. Function -> val -> SRTree val
Uni Function
f (a
t, b
2b -> b -> b
forall a. Num a => a -> a -> a
*b
iyb -> b -> b
forall a. Num a => a -> a -> a
+b
1)

    convert :: SRTree [(a, (a, Int, Int, Double))]
-> a -> [(a, (a, Int, Int, Double))]
convert (Var Int
ix) a
iy = [(a
iy, (a
0, Int
0, Int
ix, -Double
1))]
    convert (Const Double
x) a
iy = [(a
iy, (a
0, Int
2, -Int
1, Double
x))]
    convert (Param Int
ix) a
iy = [(a
iy, (a
0, Int
1, Int
ix, -Double
1))]
    convert (Uni Function
f [(a, (a, Int, Int, Double))]
t) a
iy = (a
iy, (a
1, Function -> Int
forall a. Enum a => a -> Int
fromEnum Function
f, -Int
1, -Double
1)) (a, (a, Int, Int, Double))
-> [(a, (a, Int, Int, Double))] -> [(a, (a, Int, Int, Double))]
forall a. a -> [a] -> [a]
: [(a, (a, Int, Int, Double))]
t
    convert (Bin Op
op [(a, (a, Int, Int, Double))]
l [(a, (a, Int, Int, Double))]
r) a
iy = (a
iy, (a
2, Op -> Int
forall a. Enum a => a -> Int
fromEnum Op
op, -Int
1, -Double
1)) (a, (a, Int, Int, Double))
-> [(a, (a, Int, Int, Double))] -> [(a, (a, Int, Int, Double))]
forall a. a -> [a] -> [a]
: ([(a, (a, Int, Int, Double))]
l [(a, (a, Int, Int, Double))]
-> [(a, (a, Int, Int, Double))] -> [(a, (a, Int, Int, Double))]
forall a. Semigroup a => a -> a -> a
<> [(a, (a, Int, Int, Double))]
r)

-- | minimizes the negative log-likelihood of the expression
minimizeNLL :: Distribution -> Maybe Double -> Int -> SRMatrix -> PVector -> Fix SRTree -> PVector -> (PVector, Double)
minimizeNLL :: Distribution
-> Maybe Double
-> Int
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (PVector, Double)
minimizeNLL Distribution
dist Maybe Double
msErr Int
niter SRMatrix
xss PVector
ys Fix SRTree
tree PVector
t0
  | Int
niter Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = (PVector
t0, Double
f)
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0     = (PVector
t0, Double
f)
  | Bool
otherwise  = (Comp -> Vector Double -> PVector
forall e. Comp -> Vector e -> Vector S e
fromStorableVector Comp
compMode Vector Double
t_opt, Double
f)
  where
    tree' :: Fix SRTree
tree'      = Fix SRTree -> Fix SRTree
relabelParams Fix SRTree
tree -- $ fst $ floatConstsToParam tree
    t0' :: Vector Double
t0'        = PVector -> Vector Double
forall ix e. Index ix => Array S ix e -> Vector e
toStorableVector PVector
t0
    treeArr :: [(Int, (Int, Int, Int, Double))]
treeArr    = IntMap (Int, Int, Int, Double) -> [(Int, (Int, Int, Int, Double))]
forall a. IntMap a -> [(Int, a)]
IntMap.toAscList (IntMap (Int, Int, Int, Double)
 -> [(Int, (Int, Int, Int, Double))])
-> IntMap (Int, Int, Int, Double)
-> [(Int, (Int, Int, Int, Double))]
forall a b. (a -> b) -> a -> b
$ Fix SRTree -> IntMap (Int, Int, Int, Double)
tree2arr Fix SRTree
tree'
    j2ix :: IntMap Int
j2ix       = [(Int, Int)] -> IntMap Int
forall a. [(Int, a)] -> IntMap a
IntMap.fromList ([(Int, Int)] -> IntMap Int) -> [(Int, Int)] -> IntMap Int
forall a b. (a -> b) -> a -> b
$ [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
Prelude.zip (((Int, (Int, Int, Int, Double)) -> Int)
-> [(Int, (Int, Int, Int, Double))] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
Prelude.map (Int, (Int, Int, Int, Double)) -> Int
forall a b. (a, b) -> a
fst [(Int, (Int, Int, Int, Double))]
treeArr) [Int
0..]
    (Sz Int
n)     = PVector -> Sz Int
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
size PVector
t0
    (Sz Int
m)     = PVector -> Sz Int
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
size PVector
ys
    funAndGrad :: Vector Double -> (Double, Vector Double)
funAndGrad = (Array D Int Double -> Vector Double)
-> (Double, Array D Int Double) -> (Double, Vector Double)
forall b c a. (b -> c) -> (a, b) -> (a, c)
forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second (PVector -> Vector Double
forall ix e. Index ix => Array S ix e -> Vector e
toStorableVector (PVector -> Vector Double)
-> (Array D Int Double -> PVector)
-> Array D Int Double
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. S -> Array D Int Double -> PVector
forall r e r' ix.
(Manifest r e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
computeAs S
S) ((Double, Array D Int Double) -> (Double, Vector Double))
-> (Vector Double -> (Double, Array D Int Double))
-> Vector Double
-> (Double, Vector Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> [(Int, (Int, Int, Int, Double))]
-> IntMap Int
-> PVector
-> (Double, Array D Int Double)
gradNLLArr Distribution
dist Maybe Double
msErr SRMatrix
xss PVector
ys [(Int, (Int, Int, Int, Double))]
treeArr IntMap Int
j2ix (PVector -> (Double, Array D Int Double))
-> (Vector Double -> PVector)
-> Vector Double
-> (Double, Array D Int Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Comp -> Vector Double -> PVector
forall e. Comp -> Vector e -> Vector S e
fromStorableVector Comp
compMode
    (Double
f, Array D Int Double
_)     = Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> [(Int, (Int, Int, Int, Double))]
-> IntMap Int
-> PVector
-> (Double, Array D Int Double)
gradNLLArr Distribution
dist Maybe Double
msErr SRMatrix
xss PVector
ys [(Int, (Int, Int, Int, Double))]
treeArr IntMap Int
j2ix PVector
t0 -- if there's no parameter or no iterations

    algorithm :: LocalAlgorithm
algorithm  = (Vector Double -> (Double, Vector Double))
-> Maybe VectorStorage -> LocalAlgorithm
LBFGS Vector Double -> (Double, Vector Double)
funAndGrad Maybe VectorStorage
forall a. Maybe a
Nothing
    stop :: NonEmpty StoppingCondition
stop       = Double -> StoppingCondition
ObjectiveRelativeTolerance Double
1e-6 StoppingCondition
-> [StoppingCondition] -> NonEmpty StoppingCondition
forall a. a -> [a] -> NonEmpty a
:| [Word -> StoppingCondition
MaximumEvaluations (Int -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
niter)]
    problem :: LocalProblem
problem    = Word
-> NonEmpty StoppingCondition -> LocalAlgorithm -> LocalProblem
LocalProblem (Int -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) NonEmpty StoppingCondition
stop LocalAlgorithm
algorithm
    t_opt :: Vector Double
t_opt      = case LocalProblem -> Vector Double -> Either Result Solution
minimizeLocal LocalProblem
problem Vector Double
t0' of
                  Right Solution
sol -> Solution -> Vector Double
solutionParams Solution
sol
                  Left Result
e    -> Vector Double
t0'

-- | minimizes the likelihood assuming repeating parameters in the expression 
minimizeNLLNonUnique :: Distribution -> Maybe Double -> Int -> SRMatrix -> PVector -> Fix SRTree -> PVector -> (PVector, Double)
minimizeNLLNonUnique :: Distribution
-> Maybe Double
-> Int
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (PVector, Double)
minimizeNLLNonUnique Distribution
dist Maybe Double
msErr Int
niter SRMatrix
xss PVector
ys Fix SRTree
tree PVector
t0
  | Int
niter Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = (PVector
t0, Double
f)
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0     = (PVector
t0, Double
f)
  | Bool
otherwise  = (Comp -> Vector Double -> PVector
forall e. Comp -> Vector e -> Vector S e
fromStorableVector Comp
compMode Vector Double
t_opt, Double
f)
  where
    t0' :: Vector Double
t0'        = PVector -> Vector Double
forall ix e. Index ix => Array S ix e -> Vector e
toStorableVector PVector
t0
    (Sz Int
n)     = PVector -> Sz Int
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
size PVector
t0
    (Sz Int
m)     = PVector -> Sz Int
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
size PVector
ys
    funAndGrad :: Vector Double -> (Double, Vector Double)
funAndGrad = (Array D Int Double -> Vector Double)
-> (Double, Array D Int Double) -> (Double, Vector Double)
forall b c a. (b -> c) -> (a, b) -> (a, c)
forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second (PVector -> Vector Double
forall ix e. Index ix => Array S ix e -> Vector e
toStorableVector (PVector -> Vector Double)
-> (Array D Int Double -> PVector)
-> Array D Int Double
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. S -> Array D Int Double -> PVector
forall r e r' ix.
(Manifest r e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
computeAs S
S) ((Double, Array D Int Double) -> (Double, Vector Double))
-> (Vector Double -> (Double, Array D Int Double))
-> Vector Double
-> (Double, Vector Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (Double, Array D Int Double)
gradNLLNonUnique Distribution
dist Maybe Double
msErr SRMatrix
xss PVector
ys Fix SRTree
tree (PVector -> (Double, Array D Int Double))
-> (Vector Double -> PVector)
-> Vector Double
-> (Double, Array D Int Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Comp -> Vector Double -> PVector
forall e. Comp -> Vector e -> Vector S e
fromStorableVector Comp
compMode
    (Double
f, Array D Int Double
_)     = Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (Double, Array D Int Double)
gradNLLNonUnique Distribution
dist Maybe Double
msErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
t0 -- if there's no parameter or no iterations

    algorithm :: LocalAlgorithm
algorithm  = (Vector Double -> (Double, Vector Double))
-> Maybe VectorStorage -> LocalAlgorithm
LBFGS Vector Double -> (Double, Vector Double)
funAndGrad Maybe VectorStorage
forall a. Maybe a
Nothing
    stop :: NonEmpty StoppingCondition
stop       = Double -> StoppingCondition
ObjectiveRelativeTolerance Double
1e-5 StoppingCondition
-> [StoppingCondition] -> NonEmpty StoppingCondition
forall a. a -> [a] -> NonEmpty a
:| [Word -> StoppingCondition
MaximumEvaluations (Int -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
niter)]
    problem :: LocalProblem
problem    = Word
-> NonEmpty StoppingCondition -> LocalAlgorithm -> LocalProblem
LocalProblem (Int -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) NonEmpty StoppingCondition
stop LocalAlgorithm
algorithm
    t_opt :: Vector Double
t_opt      = case LocalProblem -> Vector Double -> Either Result Solution
minimizeLocal LocalProblem
problem Vector Double
t0' of
                  Right Solution
sol -> Solution -> Vector Double
solutionParams Solution
sol
                  Left Result
e    -> Vector Double
t0'

-- | minimizes the function while keeping the parameter ix fixed (used to calculate the profile)
minimizeNLLWithFixedParam :: Distribution -> Maybe Double -> Int -> SRMatrix -> PVector -> Fix SRTree -> Int -> PVector -> PVector
minimizeNLLWithFixedParam :: Distribution
-> Maybe Double
-> Int
-> SRMatrix
-> PVector
-> Fix SRTree
-> Int
-> PVector
-> PVector
minimizeNLLWithFixedParam Distribution
dist Maybe Double
msErr Int
niter SRMatrix
xss PVector
ys Fix SRTree
tree Int
ix PVector
t0
  | Int
niter Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = PVector
t0
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0     = PVector
t0
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
m      = PVector
t0
  | Bool
otherwise  = Comp -> Vector Double -> PVector
forall e. Comp -> Vector e -> Vector S e
fromStorableVector Comp
compMode Vector Double
t_opt
  where
    t0' :: Vector Double
t0'        = PVector -> Vector Double
forall ix e. Index ix => Array S ix e -> Vector e
toStorableVector PVector
t0
    (Sz Int
n)     = PVector -> Sz Int
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
size PVector
t0
    (Sz Int
m)     = PVector -> Sz Int
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
size PVector
ys
    setTo0 :: Vector Double -> Vector Double
setTo0     = (Vector Double -> [(Int, Double)] -> Vector Double
forall a. Storable a => Vector a -> [(Int, a)] -> Vector a
VS.// [(Int
ix, Double
0.0)])
    funAndGrad :: Vector Double -> (Double, Vector Double)
funAndGrad = (Array D Int Double -> Vector Double)
-> (Double, Array D Int Double) -> (Double, Vector Double)
forall b c a. (b -> c) -> (a, b) -> (a, c)
forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second (Vector Double -> Vector Double
setTo0 (Vector Double -> Vector Double)
-> (Array D Int Double -> Vector Double)
-> Array D Int Double
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. PVector -> Vector Double
forall ix e. Index ix => Array S ix e -> Vector e
toStorableVector (PVector -> Vector Double)
-> (Array D Int Double -> PVector)
-> Array D Int Double
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. S -> Array D Int Double -> PVector
forall r e r' ix.
(Manifest r e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
computeAs S
S)((Double, Array D Int Double) -> (Double, Vector Double))
-> (Vector Double -> (Double, Array D Int Double))
-> Vector Double
-> (Double, Vector Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (Double, Array D Int Double)
gradNLLNonUnique Distribution
dist Maybe Double
msErr SRMatrix
xss PVector
ys Fix SRTree
tree (PVector -> (Double, Array D Int Double))
-> (Vector Double -> PVector)
-> Vector Double
-> (Double, Array D Int Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Comp -> Vector Double -> PVector
forall e. Comp -> Vector e -> Vector S e
fromStorableVector Comp
compMode
    (Double
f, Array D Int Double
_)     = Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (Double, Array D Int Double)
gradNLLNonUnique Distribution
dist Maybe Double
msErr SRMatrix
xss PVector
ys Fix SRTree
tree PVector
t0 -- if there's no parameter or no iterations

    algorithm :: LocalAlgorithm
algorithm  = (Vector Double -> (Double, Vector Double))
-> Maybe VectorStorage -> LocalAlgorithm
LBFGS Vector Double -> (Double, Vector Double)
funAndGrad Maybe VectorStorage
forall a. Maybe a
Nothing
    stop :: NonEmpty StoppingCondition
stop       = Double -> StoppingCondition
ObjectiveRelativeTolerance Double
1e-5 StoppingCondition
-> [StoppingCondition] -> NonEmpty StoppingCondition
forall a. a -> [a] -> NonEmpty a
:| [Word -> StoppingCondition
MaximumEvaluations (Int -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
niter)]
    problem :: LocalProblem
problem    = Word
-> NonEmpty StoppingCondition -> LocalAlgorithm -> LocalProblem
LocalProblem (Int -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) NonEmpty StoppingCondition
stop LocalAlgorithm
algorithm
    t_opt :: Vector Double
t_opt      = case LocalProblem -> Vector Double -> Either Result Solution
minimizeLocal LocalProblem
problem Vector Double
t0' of
                  Right Solution
sol -> Solution -> Vector Double
solutionParams Solution
sol
                  Left Result
e    -> Vector Double
t0'

-- | minimizes using Gaussian likelihood 
minimizeGaussian :: Int -> SRMatrix -> PVector -> Fix SRTree -> PVector -> (PVector, Double)
minimizeGaussian :: Int
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (PVector, Double)
minimizeGaussian = Distribution
-> Maybe Double
-> Int
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (PVector, Double)
minimizeNLL Distribution
Gaussian Maybe Double
forall a. Maybe a
Nothing

-- | minimizes using Binomial likelihood 
minimizeBinomial :: Int -> SRMatrix -> PVector -> Fix SRTree -> PVector -> (PVector, Double)
minimizeBinomial :: Int
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (PVector, Double)
minimizeBinomial = Distribution
-> Maybe Double
-> Int
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (PVector, Double)
minimizeNLL Distribution
Bernoulli Maybe Double
forall a. Maybe a
Nothing

-- | minimizes using Poisson likelihood 
minimizePoisson :: Int -> SRMatrix -> PVector -> Fix SRTree -> PVector -> (PVector, Double)
minimizePoisson :: Int
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (PVector, Double)
minimizePoisson = Distribution
-> Maybe Double
-> Int
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (PVector, Double)
minimizeNLL Distribution
Poisson Maybe Double
forall a. Maybe a
Nothing

-- estimates the standard error if not provided 
estimateSErr :: Distribution -> Maybe Double -> SRMatrix -> PVector -> PVector -> Fix SRTree -> Int -> Maybe Double
estimateSErr :: Distribution
-> Maybe Double
-> SRMatrix
-> PVector
-> PVector
-> Fix SRTree
-> Int
-> Maybe Double
estimateSErr Distribution
Gaussian Maybe Double
Nothing  SRMatrix
xss PVector
ys PVector
theta0 Fix SRTree
t Int
nIter = Double -> Maybe Double
forall a. a -> Maybe a
Just Double
err
  where
    theta :: PVector
theta  = (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
-> Int
-> SRMatrix
-> PVector
-> Fix SRTree
-> PVector
-> (PVector, Double)
minimizeNLL Distribution
Gaussian (Double -> Maybe Double
forall a. a -> Maybe a
Just Double
1) Int
nIter SRMatrix
xss PVector
ys Fix SRTree
t PVector
theta0
    (Sz Int
m) = PVector -> Sz Int
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
size PVector
ys
    (Sz Int
p) = PVector -> Sz Int
forall r ix e. Size r => Array r ix e -> Sz ix
forall ix e. Array S ix e -> Sz ix
size PVector
theta
    ssr :: Double
ssr    = SRMatrix -> PVector -> Fix SRTree -> PVector -> Double
sse SRMatrix
xss PVector
ys Fix SRTree
t PVector
theta
    err :: Double
err    = 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
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
p)
estimateSErr Distribution
_        (Just Double
s) SRMatrix
_   PVector
_  PVector
_ Fix SRTree
_ Int
_   = Double -> Maybe Double
forall a. a -> Maybe a
Just Double
s
estimateSErr Distribution
_        Maybe Double
_        SRMatrix
_   PVector
_  PVector
_ Fix SRTree
_ Int
_   = Maybe Double
forall a. Maybe a
Nothing