{-|
Module      : MachineLearning.Model.Regression
Description : TIR expression data structures
Copyright   : (c) Fabricio Olivetti de Franca, 2022
License     : GPL-3
Maintainer  : fabricio.olivetti@gmail.com
Stability   : experimental
Portability : POSIX

Fitting functions for the coefficients.
-}
module MachineLearning.Model.Regression 
  ( fitTask
  , predictTask
  , evalPenalty
  , applyMeasures
  , nonlinearFit
  , tirToMatrix
  ) where

import Data.Bifunctor
import           Data.SRTree                   (evalFun, Function(..), OptIntPow(..), evalFun, inverseFunc)
import MachineLearning.TIR       (TIR(..),  Individual(..), Dataset, Constraint, assembleTree, replaceConsts)
import qualified Data.Vector           as V
import qualified Data.Vector.Storable  as VS
import Data.Vector           ((!))
import           Data.List                     (nub)
import           Data.Vector.Storable          (Vector, splitAt)
import           Numeric.LinearAlgebra         ((<\>), Matrix)
import           Numeric.GSL.Fitting           (nlFitting, FittingMethod(..))
import           Prelude               hiding  (splitAt)

import qualified Numeric.LinearAlgebra                     as LA

import MachineLearning.TIR (Individual(..))
import MachineLearning.Utils.Config (Task(..), Penalty(..))
import MachineLearning.Model.Measure (Measure(..))

-- * IT specific stuff

createQ :: Vector Double -> LA.Matrix Double -> LA.Matrix Double
createQ :: Vector Double -> Matrix Double -> Matrix Double
createQ Vector Double
ys = [Vector Double] -> Matrix Double
forall t. Element t => [Vector t] -> Matrix t
LA.fromColumns ([Vector Double] -> Matrix Double)
-> (Matrix Double -> [Vector Double])
-> Matrix Double
-> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector Double -> Vector Double)
-> [Vector Double] -> [Vector Double]
forall a b. (a -> b) -> [a] -> [b]
map (Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
*Vector Double
ys) ([Vector Double] -> [Vector Double])
-> (Matrix Double -> [Vector Double])
-> Matrix Double
-> [Vector Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix Double -> [Vector Double]
forall t. Element t => Matrix t -> [Vector t]
LA.toColumns
{-# INLINE createQ #-}

avg :: Vector Double -> Vector Double
avg :: Vector Double -> Vector Double
avg Vector Double
ys = [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
LA.fromList [Vector Double -> Double
forall (c :: * -> *) e. Container c e => c e -> e
LA.sumElements Vector Double
ys Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Vector Double -> IndexOf Vector
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
LA.size Vector Double
ys)]
{-# INLINE avg #-}

-- | transform a data matrix using a TIR expression. This function returns
-- a tuple with the transformed data of the numerator and denominator, respectivelly.
-- Each column of the transformed data represents one term of the TIR expression.
tirToMatrix :: Dataset Double -> TIR -> (LA.Matrix Double, LA.Matrix Double)
tirToMatrix :: Dataset Double -> TIR -> (Matrix Double, Matrix Double)
tirToMatrix Dataset Double
xss (TIR Function
_ Sigma
p Sigma
q) = ([Vector Double] -> Matrix Double)
-> ([Vector Double] -> Matrix Double)
-> ([Vector Double], [Vector Double])
-> (Matrix Double, Matrix Double)
forall (p :: * -> * -> *) a b c d.
Bifunctor p =>
(a -> b) -> (c -> d) -> p a c -> p b d
bimap ([Vector Double] -> Matrix Double
forall t. Element t => [Vector t] -> Matrix t
LA.fromColumns ([Vector Double] -> Matrix Double)
-> ([Vector Double] -> [Vector Double])
-> [Vector Double]
-> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector Double
biasVector Double -> [Vector Double] -> [Vector Double]
forall a. a -> [a] -> [a]
:)) [Vector Double] -> Matrix Double
forall t. Element t => [Vector t] -> Matrix t
LA.fromColumns ([Vector Double]
p', [Vector Double]
q')
  where
    bias :: Vector Double
bias      = Dataset Double -> Vector Double
forall a. Vector a -> a
V.head Dataset Double
xss
    xss' :: Dataset Double
xss'      = Dataset Double -> Dataset Double
forall a. Vector a -> Vector a
V.tail Dataset Double
xss
    sigma2mtx :: [(a, Function, [(Int, Int)])] -> [Vector Double]
sigma2mtx = ((a, Function, [(Int, Int)]) -> Vector Double)
-> [(a, Function, [(Int, Int)])] -> [Vector Double]
forall a b. (a -> b) -> [a] -> [b]
map (\(a
_, Function
g, [(Int, Int)]
ps) -> Function -> Vector Double -> Vector Double
forall val. Floating val => Function -> val -> val
evalFun Function
g (Vector Double -> Vector Double) -> Vector Double -> Vector Double
forall a b. (a -> b) -> a -> b
$ [(Int, Int)] -> Vector Double
evalPi [(Int, Int)]
ps)
    evalPi :: [(Int, Int)] -> Vector Double
evalPi    = ((Int, Int) -> Vector Double -> Vector Double)
-> Vector Double -> [(Int, Int)] -> Vector Double
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\(Int
ix, Int
k) Vector Double
acc -> Vector Double
acc Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
* (Dataset Double
xss' Dataset Double -> Int -> Vector Double
forall a. Vector a -> Int -> a
! Int
ix Vector Double -> Int -> Vector Double
forall a b. (Fractional a, Integral b) => a -> b -> a
^^ Int
k)) Vector Double
1
    p' :: [Vector Double]
p'        = Sigma -> [Vector Double]
forall {a}. [(a, Function, [(Int, Int)])] -> [Vector Double]
sigma2mtx Sigma
p
    q' :: [Vector Double]
q'        = Sigma -> [Vector Double]
forall {a}. [(a, Function, [(Int, Int)])] -> [Vector Double]
sigma2mtx Sigma
q
{-# INLINE tirToMatrix #-}

-- | Fits a linear model using l2-penalty
ridge :: Matrix Double -> Vector Double -> Matrix Double
ridge :: Matrix Double -> Vector Double -> Matrix Double
ridge Matrix Double
a Vector Double
b = Matrix Double
oA Matrix Double -> Matrix Double -> Matrix Double
forall (c :: * -> *) t.
(LSDiv c, Field t) =>
Matrix t -> c t -> c t
<\> Matrix Double
oB
  where
   mu :: Matrix Double
mu = Matrix Double
0.01

   a' :: Matrix Double
a' = Matrix Double -> Matrix Double
forall m mt. Transposable m mt => m -> mt
LA.tr Matrix Double
a
   b' :: Matrix Double
b' = Matrix Double -> Matrix Double
forall m mt. Transposable m mt => m -> mt
LA.tr (Matrix Double -> Matrix Double) -> Matrix Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Vector Double -> Matrix Double
forall a. Storable a => Vector a -> Matrix a
LA.asColumn Vector Double
b
   oA :: Matrix Double
oA = (Matrix Double
a' Matrix Double -> Matrix Double -> Matrix Double
forall a. Semigroup a => a -> a -> a
<> Matrix Double -> Matrix Double
forall m mt. Transposable m mt => m -> mt
LA.tr Matrix Double
a') Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
+ (Matrix Double
mu Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
* Int -> Matrix Double
forall a. (Num a, Element a) => Int -> Matrix a
LA.ident (Matrix Double -> Int
forall t. Matrix t -> Int
LA.rows Matrix Double
a'))
   oB :: Matrix Double
oB = Matrix Double
a' Matrix Double -> Matrix Double -> Matrix Double
forall a. Semigroup a => a -> a -> a
<> Matrix Double -> Matrix Double
forall m mt. Transposable m mt => m -> mt
LA.tr Matrix Double
b'
{-# INLINE ridge #-}

-- | Predicts a linear model
predict :: Matrix Double -> Vector Double -> Vector Double 
predict :: Matrix Double -> Vector Double -> Vector Double
predict Matrix Double
xs Vector Double
w | Matrix Double -> Int
forall t. Matrix t -> Int
LA.cols Matrix Double
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Vector Double -> IndexOf Vector
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
LA.size Vector Double
w = Matrix Double
xs Matrix Double -> Vector Double -> Vector Double
forall t. Numeric t => Matrix t -> Vector t -> Vector t
LA.#> Vector Double
w
             | Bool
otherwise = [Char] -> Vector Double
forall a. HasCallStack => [Char] -> a
error ([Char] -> Vector Double) -> [Char] -> Vector Double
forall a b. (a -> b) -> a -> b
$ [Char]
"predict: " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ (Int, Int) -> [Char]
forall a. Show a => a -> [Char]
show (Matrix Double -> IndexOf Matrix
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
LA.size Matrix Double
xs) [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show (Vector Double -> IndexOf Vector
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
LA.size Vector Double
w)
{-# INLINE predict #-}

-- | Solve the OLS *zss*w = ys*
solveOLS :: Matrix Double -> Vector Double -> Vector Double 
solveOLS :: Matrix Double -> Vector Double -> Vector Double
solveOLS Matrix Double
zss Vector Double
ys = Matrix Double
zss Matrix Double -> Vector Double -> Vector Double
forall (c :: * -> *) t.
(LSDiv c, Field t) =>
Matrix t -> c t -> c t
<\> Vector Double
ys
{-# INLINE solveOLS #-}

-- | Applies OLS and returns a Solution
-- if the expression is invalid, it returns Infinity as a fitness
--regress :: Matrix Double -> Vector Double -> [Vector Double]
regress :: TIR -> Dataset Double -> Vector Double -> [Vector Double]
regress :: TIR -> Dataset Double -> Vector Double -> [Vector Double]
regress TIR
tir Dataset Double
xss Vector Double
ys = [Vector Double
ws]
  where
    (Matrix Double
zssP, Matrix Double
zssQ) = Dataset Double -> TIR -> (Matrix Double, Matrix Double)
tirToMatrix Dataset Double
xss TIR
tir
    ys' :: Vector Double
ys'          = Function -> Vector Double -> Vector Double
forall val. Floating val => Function -> val -> val
evalFun (Function -> Function
inverseFunc (Function -> Function) -> Function -> Function
forall a b. (a -> b) -> a -> b
$ TIR -> Function
_funY TIR
tir) Vector Double
ys
    zssQy :: Matrix Double
zssQy        = Vector Double -> Matrix Double -> Matrix Double
createQ Vector Double
ys' Matrix Double
zssQ
    zss :: Matrix Double
zss          = if Matrix Double -> Int
forall t. Matrix t -> Int
LA.cols Matrix Double
zssQ Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
1
                       then Matrix Double
zssP Matrix Double -> Matrix Double -> Matrix Double
forall t. Element t => Matrix t -> Matrix t -> Matrix t
LA.||| Matrix Double -> Matrix Double
forall a. Num a => a -> a
negate Matrix Double
zssQy 
                       else Matrix Double
zssP
    ws :: Vector Double
ws             = if Matrix Double -> Int
forall t. Matrix t -> Int
LA.cols Matrix Double
zss Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1
                       then Vector Double -> Vector Double
avg Vector Double
ys'
                       else Matrix Double -> Vector Double -> Vector Double
solveOLS Matrix Double
zss Vector Double
ys'
-- regress zss ys = [solveOLS zss ys]
{-# INLINE regress #-}

-- | Applies conjugate gradient for binary classification
--classify :: Matrix Double -> Vector Double -> [Vector Double]
classify :: Int -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
classify :: Int -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
classify Int
niter TIR
tir Dataset Double
xss Vector Double
ys = [Vector Double
ws]
  where
    ws :: Vector Double
ws           = Int
-> Matrix Double
-> Matrix Double
-> Vector Double
-> (Vector Double -> Vector Double)
-> (Vector Double -> Vector Double)
-> Vector Double
-> Vector Double
nonlinearFit Int
niter Matrix Double
zssP Matrix Double
zssQ Vector Double
ys Vector Double -> Vector Double
forall a. Floating a => a -> a
sigmoid Vector Double -> Vector Double
forall a. Floating a => a -> a
dsigmoid Vector Double
theta0
    theta0 :: Vector Double
theta0       = Double -> Int -> Vector Double
forall e d (c :: * -> *). Konst e d c => e -> d -> c e
LA.konst Double
0 (Matrix Double -> Int
forall t. Matrix t -> Int
LA.cols Matrix Double
zssP Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Matrix Double -> Int
forall t. Matrix t -> Int
LA.cols Matrix Double
zssQ)
    (Matrix Double
zssP, Matrix Double
zssQ) = Dataset Double -> TIR -> (Matrix Double, Matrix Double)
tirToMatrix Dataset Double
xss TIR
tir

-- | Applies conjugate gradient for one-vs-all classification
--classifyMult :: Matrix Double -> Vector Double -> [Vector Double]
classifyMult :: Int -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
classifyMult :: Int -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
classifyMult Int
niter TIR
tir Dataset Double
xss Vector Double
ys = (Vector Double -> Vector Double -> Vector Double)
-> [Vector Double] -> [Vector Double] -> [Vector Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Vector Double -> Vector Double -> Vector Double
minimize [Vector Double]
yss [Vector Double]
theta0
  where
    numLabels :: Int
numLabels    = [Double] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Double] -> Int) -> [Double] -> Int
forall a b. (a -> b) -> a -> b
$ [Double] -> [Double]
forall a. Eq a => [a] -> [a]
nub ([Double] -> [Double]) -> [Double] -> [Double]
forall a b. (a -> b) -> a -> b
$ Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
LA.toList Vector Double
ys
    yss :: [Vector Double]
yss          = (Int -> Vector Double) -> [Int] -> [Vector Double]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Vector Double
forall {b} {a}. (Storable b, Integral a, Num b) => a -> Vector b
f [Int
0 .. Int
numLabelsInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]
    minimize :: Vector Double -> Vector Double -> Vector Double
minimize Vector Double
y Vector Double
t = Int
-> Matrix Double
-> Matrix Double
-> Vector Double
-> (Vector Double -> Vector Double)
-> (Vector Double -> Vector Double)
-> Vector Double
-> Vector Double
nonlinearFit Int
niter Matrix Double
zssP Matrix Double
zssQ Vector Double
y Vector Double -> Vector Double
forall a. Floating a => a -> a
sigmoid Vector Double -> Vector Double
forall a. Floating a => a -> a
dsigmoid Vector Double
t
    theta0 :: [Vector Double]
theta0       = Int -> Vector Double -> [Vector Double]
forall a. Int -> a -> [a]
replicate Int
numLabels (Vector Double -> [Vector Double])
-> Vector Double -> [Vector Double]
forall a b. (a -> b) -> a -> b
$ Double -> Int -> Vector Double
forall e d (c :: * -> *). Konst e d c => e -> d -> c e
LA.konst Double
0 (Matrix Double -> Int
forall t. Matrix t -> Int
LA.cols Matrix Double
zssP Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Matrix Double -> Int
forall t. Matrix t -> Int
LA.cols Matrix Double
zssQ)
    (Matrix Double
zssP, Matrix Double
zssQ) = Dataset Double -> TIR -> (Matrix Double, Matrix Double)
tirToMatrix Dataset Double
xss TIR
tir
    
    f :: a -> Vector b
f a
sample = (Double -> b) -> Vector Double -> Vector b
forall a b.
(Storable a, Storable b) =>
(a -> b) -> Vector a -> Vector b
VS.map (\Double
a -> if Double -> a
forall a b. (RealFrac a, Integral b) => a -> b
round Double
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
sample then b
1 else b
0) Vector Double
ys

-- | chooses the appropriate fitting function
--fitTask :: Task -> Matrix Double -> Vector Double -> [Vector Double]
fitTask :: Task -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
fitTask :: Task -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
fitTask Task
Regression             = TIR -> Dataset Double -> Vector Double -> [Vector Double]
regress
fitTask (RegressionNL Int
niter)   = Int -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
regressNL Int
niter
fitTask (Classification Int
niter) = Int -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
classify Int
niter
fitTask (ClassMult Int
niter)      = Int -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
classifyMult Int
niter
{-# INLINE fitTask #-}

-- | sigmoid function for classification.
sigmoid :: Floating a => a -> a
sigmoid :: forall a. Floating a => a -> a
sigmoid a
z = a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
1a -> a -> a
forall a. Num a => a -> a -> a
+a -> a
forall a. Floating a => a -> a
exp(-a
z))
{-# INLINE sigmoid #-}

-- | derivative sigmoid function for classification.
dsigmoid :: Floating a => a -> a
dsigmoid :: forall a. Floating a => a -> a
dsigmoid a
z = a -> a
forall a. Floating a => a -> a
sigmoid a
z a -> a -> a
forall a. Num a => a -> a -> a
* (a
1 a -> a -> a
forall a. Num a => a -> a -> a
- a -> a
forall a. Floating a => a -> a
sigmoid a
z)
{-# INLINE dsigmoid #-}

-- | chooses the appropriate prediction function
predictTask :: Task -> [Vector Double] -> Vector Double
predictTask :: Task -> [Vector Double] -> Vector Double
predictTask Task
_ []                   = [Char] -> Vector Double
forall a. HasCallStack => [Char] -> a
error [Char]
"predictTask: empty coefficients matrix"
predictTask Task
Regression [Vector Double]
yss         = [Vector Double] -> Vector Double
forall a. [a] -> a
head [Vector Double]
yss
predictTask (RegressionNL Int
_) [Vector Double]
yss   = [Vector Double] -> Vector Double
forall a. [a] -> a
head [Vector Double]
yss
predictTask (Classification Int
_) [Vector Double]
yss = Vector Double -> Vector Double
forall a. Floating a => a -> a
sigmoid (Vector Double -> Vector Double) -> Vector Double -> Vector Double
forall a b. (a -> b) -> a -> b
$ [Vector Double] -> Vector Double
forall a. [a] -> a
head [Vector Double]
yss
predictTask (ClassMult Int
_) [Vector Double]
yss      = [Double] -> Vector Double
LA.vector ([Double] -> Vector Double) -> [Double] -> Vector Double
forall a b. (a -> b) -> a -> b
$ (Vector Double -> Double) -> [Vector Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Double)
-> (Vector Double -> Int) -> Vector Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> Int
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
LA.maxIndex) ([Vector Double] -> [Double]) -> [Vector Double] -> [Double]
forall a b. (a -> b) -> a -> b
$ Matrix Double -> [Vector Double]
forall t. Element t => Matrix t -> [Vector t]
LA.toRows (Matrix Double -> [Vector Double])
-> Matrix Double -> [Vector Double]
forall a b. (a -> b) -> a -> b
$ [Vector Double] -> Matrix Double
forall t. Element t => [Vector t] -> Matrix t
LA.fromColumns ([Vector Double] -> Matrix Double)
-> [Vector Double] -> Matrix Double
forall a b. (a -> b) -> a -> b
$ (Vector Double -> Vector Double)
-> [Vector Double] -> [Vector Double]
forall a b. (a -> b) -> [a] -> [b]
map Vector Double -> Vector Double
forall a. Floating a => a -> a
sigmoid [Vector Double]
yss
{-# INLINE predictTask #-}

-- | evals the penalty function
evalPenalty :: Penalty -> Int -> Double -> Double
evalPenalty :: Penalty -> Int -> Double -> Double
evalPenalty Penalty
NoPenalty Int
_   Double
_   = Double
0.0
evalPenalty (Len Double
c)   Int
len Double
_   = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
len Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c
evalPenalty (Shape Double
c) Int
_   Double
val = Double
valDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
c
{-# INLINE evalPenalty #-}

-- | applies a list of performance measures
applyMeasures :: [Measure] -> Vector Double -> Vector Double -> [Double]
applyMeasures :: [Measure] -> Vector Double -> Vector Double -> [Double]
applyMeasures [Measure]
measures Vector Double
ys Vector Double
ysHat = (Measure -> Double) -> [Measure] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (((Vector Double -> Vector Double -> Double)
-> (Vector Double, Vector Double) -> Double
forall a b c. (a -> b -> c) -> (a, b) -> c
`uncurry` (Vector Double
ys, Vector Double
ysHat)) ((Vector Double -> Vector Double -> Double) -> Double)
-> (Measure -> Vector Double -> Vector Double -> Double)
-> Measure
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Measure -> Vector Double -> Vector Double -> Double
_fun) [Measure]
measures
{-# INLINE applyMeasures #-}

regressNL :: Int -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
regressNL :: Int -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
regressNL Int
niter TIR
tir Dataset Double
xss Vector Double
ys = [Vector Double
ws]
  where
    ws :: Vector Double
ws           = Int
-> Matrix Double
-> Matrix Double
-> Vector Double
-> (Vector Double -> Vector Double)
-> (Vector Double -> Vector Double)
-> Vector Double
-> Vector Double
nonlinearFit Int
niter Matrix Double
zssP Matrix Double
zssQ Vector Double
ys Vector Double -> Vector Double
f Vector Double -> Vector Double
f' Vector Double
theta0
    f :: Vector Double -> Vector Double
f            = Function -> Vector Double -> Vector Double
forall val. Floating val => Function -> val -> val
evalFun (Function -> Vector Double -> Vector Double)
-> Function -> Vector Double -> Vector Double
forall a b. (a -> b) -> a -> b
$ TIR -> Function
_funY TIR
tir
    f' :: Vector Double -> Vector Double
f'           = Function -> Vector Double -> Vector Double
forall val. (Eq val, Floating val) => Function -> val -> val
derivative (Function -> Vector Double -> Vector Double)
-> Function -> Vector Double -> Vector Double
forall a b. (a -> b) -> a -> b
$ TIR -> Function
_funY TIR
tir
    theta0 :: Vector Double
theta0       = [Vector Double] -> Vector Double
forall a. [a] -> a
head ([Vector Double] -> Vector Double)
-> [Vector Double] -> Vector Double
forall a b. (a -> b) -> a -> b
$ TIR -> Dataset Double -> Vector Double -> [Vector Double]
regress TIR
tir Dataset Double
xss Vector Double
ys
    (Matrix Double
zssP, Matrix Double
zssQ) = Dataset Double -> TIR -> (Matrix Double, Matrix Double)
tirToMatrix Dataset Double
xss TIR
tir

-- | Non-linear optimization using Levenberg-Marquardt method.
--nonlinearFit :: Monad m => Vector Double -> Matrix Double -> Matrix Double -> Vector Double -> m (Vector Double)
nonlinearFit :: Int
             -> Matrix Double 
             -> Matrix Double 
             -> Vector Double 
             -> (Vector Double -> Vector Double) 
             -> (Vector Double -> Vector Double) 
             -> Vector Double 
             -> Vector Double
nonlinearFit :: Int
-> Matrix Double
-> Matrix Double
-> Vector Double
-> (Vector Double -> Vector Double)
-> (Vector Double -> Vector Double)
-> Vector Double
-> Vector Double
nonlinearFit Int
niter Matrix Double
zssP Matrix Double
zssQ Vector Double
ys Vector Double -> Vector Double
f Vector Double -> Vector Double
f' Vector Double
theta0 = (Vector Double, Matrix Double) -> Vector Double
forall a b. (a, b) -> a
fst ((Vector Double, Matrix Double) -> Vector Double)
-> (Vector Double, Matrix Double) -> Vector Double
forall a b. (a -> b) -> a -> b
$ FittingMethod
-> Double
-> Double
-> Int
-> (Vector Double -> Vector Double)
-> (Vector Double -> Matrix Double)
-> Vector Double
-> (Vector Double, Matrix Double)
nlFitting FittingMethod
LevenbergMarquardtScaled Double
1e-6 Double
1e-6 Int
niter Vector Double -> Vector Double
model' Vector Double -> Matrix Double
jacob' Vector Double
theta0
  where
    model' :: Vector Double -> Vector Double
model'       = (Vector Double -> Vector Double)
-> Vector Double
-> Matrix Double
-> Matrix Double
-> Vector Double
-> Vector Double
model Vector Double -> Vector Double
f Vector Double
ys Matrix Double
zssP Matrix Double
zssQ
    jacob' :: Vector Double -> Matrix Double
jacob'       = (Vector Double -> Vector Double)
-> Matrix Double -> Matrix Double -> Vector Double -> Matrix Double
jacob Vector Double -> Vector Double
f' Matrix Double
zssP Matrix Double
zssQ    

-- | calculates the error given the parameter vector beta
model :: (Vector Double -> Vector Double) -> Vector Double -> Matrix Double -> Matrix Double -> Vector Double -> Vector Double
model :: (Vector Double -> Vector Double)
-> Vector Double
-> Matrix Double
-> Matrix Double
-> Vector Double
-> Vector Double
model Vector Double -> Vector Double
f Vector Double
ys Matrix Double
zssP Matrix Double
zssQ Vector Double
beta 
  | Matrix Double -> Int
forall t. Matrix t -> Int
LA.cols Matrix Double
zssQ Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = Vector Double -> Vector Double
f Vector Double
ysHat_P Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
- Vector Double
ys
  | Bool
otherwise         = Vector Double -> Vector Double
f Vector Double
ysHat Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
- Vector Double
ys
  where
    (Vector Double
betaP, Vector Double
betaQ) = Int -> Vector Double -> (Vector Double, Vector Double)
forall a. Storable a => Int -> Vector a -> (Vector a, Vector a)
splitAt (Matrix Double -> Int
forall t. Matrix t -> Int
LA.cols Matrix Double
zssP) Vector Double
beta
    ysHat_P :: Vector Double
ysHat_P        = Matrix Double -> Vector Double -> Vector Double
predict Matrix Double
zssP Vector Double
betaP
    ysHat_Q :: Vector Double
ysHat_Q        = if Matrix Double -> Int
forall t. Matrix t -> Int
LA.cols Matrix Double
zssQ Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Vector Double
0 else Matrix Double -> Vector Double -> Vector Double
predict Matrix Double
zssQ Vector Double
betaQ
    ysHat :: Vector Double
ysHat          = Vector Double
ysHat_P Vector Double -> Vector Double -> Vector Double
forall a. Fractional a => a -> a -> a
/ (Vector Double
1 Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
+ Vector Double
ysHat_Q)

-- | calculates the Jacobian given the parameter vector beta. Doesn't support
-- the outer transformation function.
jacob :: (Vector Double -> Vector Double) -> Matrix Double -> Matrix Double -> Vector Double -> Matrix Double
jacob :: (Vector Double -> Vector Double)
-> Matrix Double -> Matrix Double -> Vector Double -> Matrix Double
jacob Vector Double -> Vector Double
f Matrix Double
zssP Matrix Double
zssQ Vector Double
beta | Matrix Double -> Int
forall t. Matrix t -> Int
LA.cols Matrix Double
zssQ Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = Matrix Double
zssP
                       | Bool
otherwise         = [Vector Double] -> Matrix Double
forall t. Element t => [Vector t] -> Matrix t
LA.fromColumns ([Vector Double] -> Matrix Double)
-> [Vector Double] -> Matrix Double
forall a b. (a -> b) -> a -> b
$ [Vector Double]
pjac [Vector Double] -> [Vector Double] -> [Vector Double]
forall a. Semigroup a => a -> a -> a
<> [Vector Double]
qjac
  where
    (Vector Double
betaP, Vector Double
betaQ) = Int -> Vector Double -> (Vector Double, Vector Double)
forall a. Storable a => Int -> Vector a -> (Vector a, Vector a)
splitAt (Matrix Double -> Int
forall t. Matrix t -> Int
LA.cols Matrix Double
zssP) Vector Double
beta
    ysHat_P :: Vector Double
ysHat_P        = Matrix Double -> Vector Double -> Vector Double
predict Matrix Double
zssP Vector Double
betaP
    ysHat_Q :: Vector Double
ysHat_Q        = Matrix Double -> Vector Double -> Vector Double
predict Matrix Double
zssQ Vector Double
betaQ
    ysHat :: Vector Double
ysHat          = Vector Double -> Vector Double
f (Vector Double -> Vector Double) -> Vector Double -> Vector Double
forall a b. (a -> b) -> a -> b
$ Vector Double
ysHat_P Vector Double -> Vector Double -> Vector Double
forall a. Fractional a => a -> a -> a
/ (Vector Double
1 Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
+ Vector Double
ysHat_Q)
    pjac :: [Vector Double]
pjac           = [Vector Double
ysHat Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
* Vector Double
c Vector Double -> Vector Double -> Vector Double
forall a. Fractional a => a -> a -> a
/ Vector Double
ysHat_Q | Vector Double
c <- Matrix Double -> [Vector Double]
forall t. Element t => Matrix t -> [Vector t]
LA.toColumns Matrix Double
zssP]
    qjac :: [Vector Double]
qjac           = [-Vector Double
ysHat Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
* Vector Double
c Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
* (Vector Double
ysHat_P Vector Double -> Vector Double -> Vector Double
forall a. Fractional a => a -> a -> a
/ (Vector Double
1 Vector Double -> Vector Double -> Vector Double
forall a. Num a => a -> a -> a
+ Vector Double
ysHat_Q)Vector Double -> Integer -> Vector Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2) | Vector Double
c <- Matrix Double -> [Vector Double]
forall t. Element t => Matrix t -> [Vector t]
LA.toColumns Matrix Double
zssQ]

derivative :: (Eq val, Floating val) => Function -> val -> val
derivative :: forall val. (Eq val, Floating val) => Function -> val -> val
derivative Function
Id      = val -> val -> val
forall a b. a -> b -> a
const val
1
derivative Function
Abs     = \val
x -> val
x val -> val -> val
forall a. Fractional a => a -> a -> a
/ val -> val
forall a. Num a => a -> a
abs val
x
derivative Function
Sin     = val -> val
forall a. Floating a => a -> a
cos
derivative Function
Cos     = val -> val
forall a. Num a => a -> a
negate(val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
.val -> val
forall a. Floating a => a -> a
sin
derivative Function
Tan     = val -> val
forall a. Fractional a => a -> a
recip (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (val -> val -> val
forall a. Floating a => a -> a -> a
**val
2.0) (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. val -> val
forall a. Floating a => a -> a
cos
derivative Function
Sinh    = val -> val
forall a. Floating a => a -> a
cosh
derivative Function
Cosh    = val -> val
forall a. Floating a => a -> a
sinh
derivative Function
Tanh    = (val
1val -> val -> val
forall a. Num a => a -> a -> a
-) (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (val -> val -> val
forall a. Floating a => a -> a -> a
**val
2.0) (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. val -> val
forall a. Floating a => a -> a
tanh
derivative Function
ASin    = val -> val
forall a. Fractional a => a -> a
recip (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. val -> val
forall a. Floating a => a -> a
sqrt (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (val
1val -> val -> val
forall a. Num a => a -> a -> a
-) (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (val -> Integer -> val
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)
derivative Function
ACos    = val -> val
forall a. Num a => a -> a
negate (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. val -> val
forall a. Fractional a => a -> a
recip (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. val -> val
forall a. Floating a => a -> a
sqrt (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (val
1val -> val -> val
forall a. Num a => a -> a -> a
-) (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (val -> Integer -> val
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)
derivative Function
ATan    = val -> val
forall a. Fractional a => a -> a
recip (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (val
1val -> val -> val
forall a. Num a => a -> a -> a
+) (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (val -> Integer -> val
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)
derivative Function
ASinh   = val -> val
forall a. Fractional a => a -> a
recip (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. val -> val
forall a. Floating a => a -> a
sqrt (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (val
1val -> val -> val
forall a. Num a => a -> a -> a
+) (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (val -> Integer -> val
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)
derivative Function
ACosh   = \val
x -> val
1 val -> val -> val
forall a. Fractional a => a -> a -> a
/ (val -> val
forall a. Floating a => a -> a
sqrt (val
xval -> val -> val
forall a. Num a => a -> a -> a
-val
1) val -> val -> val
forall a. Num a => a -> a -> a
* val -> val
forall a. Floating a => a -> a
sqrt (val
xval -> val -> val
forall a. Num a => a -> a -> a
+val
1))
derivative Function
ATanh   = val -> val
forall a. Fractional a => a -> a
recip (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (val
1val -> val -> val
forall a. Num a => a -> a -> a
-) (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (val -> Integer -> val
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)
derivative Function
Sqrt    = val -> val
forall a. Fractional a => a -> a
recip (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (val
2val -> val -> val
forall a. Num a => a -> a -> a
*) (val -> val) -> (val -> val) -> val -> val
forall b c a. (b -> c) -> (a -> b) -> a -> c
. val -> val
forall a. Floating a => a -> a
sqrt
derivative Function
Square  = (val
2val -> val -> val
forall a. Num a => a -> a -> a
*)
derivative Function
Exp     = val -> val
forall a. Floating a => a -> a
exp
derivative Function
Log     = val -> val
forall a. Fractional a => a -> a
recip
{-# INLINE derivative #-}
-- (w1 * p1 + w2 * p2) / (1 + w3 * p3 + w4 * p4)
-- d/dw1 = p1 / (1 + w3 * p3 + w4 * p4)
-- d/dw2 = p2 / (1 + w3 * p3 + w4 * p4)
-- d/dw3 = -p3 * (w1 * p2 + w2 * p2) / (1 + w3 * p3 + w4 * p4)^2
{-
toEv :: SRTree Int Double -> (V.Vector Double -> Double)
toEv (Var !ix) = (`V.unsafeIndex` ix)
toEv (Const !val) = const val
toEv (Add !l !r) = jn (+) (toEv l) (toEv r)
toEv (Mul !l !r) = jn (*) (toEv l) (toEv r)
toEv (Fun Exp !t) = exp . toEv t
{-# INLINE toEv #-}

jn :: (Double -> Double -> Double) -> (V.Vector Double -> Double) -> (V.Vector Double -> Double) -> (V.Vector Double -> Double)
jn op f g = \x -> op (f x) (g x)
{-# INLINE jn #-}
-}