{-# LANGUAGE FlexibleContexts #-}
{-|
Module      : IT.Eval
Description : Evaluation function for IT expressions.
Copyright   : (c) Fabricio Olivetti de Franca, 2020
License     : GPL-3
Maintainer  : fabricio.olivetti@gmail.com
Stability   : experimental
Portability : POSIX

Definition of functions pertaining to the conversion
and evaluation of an IT-expression. 

TODO: move interval evaluation to IT.Eval.Interval 
-}
{-# OPTIONS_GHC -Wno-unused-matches #-}

module IT.Eval where

import qualified Data.Vector as V
import qualified Data.Vector.Storable as VV
import qualified Numeric.LinearAlgebra as LA
import qualified Data.IntMap.Strict as M
import Numeric.Interval hiding (null)
import Data.Bifunctor

import IT
import IT.Algorithms

-- | log(x+1)
log1p :: Floating a => a -> a
log1p :: a -> a
log1p a
x = a -> a
forall a. Floating a => a -> a
log (a
xa -> a -> a
forall a. Num a => a -> a -> a
+a
1)

-- * Evaluation of the Transformation functions 
-- It supports any type that is an instance of Floating.

-- | Evaluate the transformation function 'f' in a data point 'x'.
transform :: Floating a => Transformation -> a -> a
transform :: Transformation -> a -> a
transform Transformation
Id      = a -> a
forall a. a -> a
id
transform Transformation
Sin     = a -> a
forall a. Floating a => a -> a
sin
transform Transformation
Cos     = a -> a
forall a. Floating a => a -> a
cos
transform Transformation
Tan     = a -> a
forall a. Floating a => a -> a
tan
transform Transformation
Tanh    = a -> a
forall a. Floating a => a -> a
tanh
transform Transformation
Sqrt    = a -> a
forall a. Floating a => a -> a
sqrt
transform Transformation
SqrtAbs = a -> a
forall a. Floating a => a -> a
sqrt(a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.a -> a
forall a. Num a => a -> a
abs
transform Transformation
Exp     = a -> a
forall a. Floating a => a -> a
exp
transform Transformation
Log     = a -> a
forall a. Floating a => a -> a
log
transform Transformation
Log1p   = a -> a
forall a. Floating a => a -> a
log1p

-- | Evaluate the derivative of a transformation function. 
derivative :: Floating a => Transformation -> a -> a
derivative :: Transformation -> a -> a
derivative Transformation
Id      = a -> a -> a
forall a b. a -> b -> a
const a
1
derivative Transformation
Sin     = a -> a
forall a. Floating a => a -> a
cos
derivative Transformation
Cos     = a -> a
forall a. Num a => a -> a
negate(a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.a -> a
forall a. Floating a => a -> a
sin
derivative Transformation
Tan     = a -> a
forall a. Fractional a => a -> a
recip (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a
forall a. Floating a => a -> a -> a
**a
2.0) (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> a
forall a. Floating a => a -> a
cos
derivative Transformation
Tanh    = (a
1a -> a -> a
forall a. Num a => a -> a -> a
-) (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a
forall a. Floating a => a -> a -> a
**a
2.0) (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> a
forall a. Floating a => a -> a
tanh
derivative Transformation
Sqrt    = a -> a
forall a. Fractional a => a -> a
recip (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a
2a -> a -> a
forall a. Num a => a -> a -> a
*) (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> a
forall a. Floating a => a -> a
sqrt
derivative Transformation
SqrtAbs = \a
x -> a
x a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
2a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Num a => a -> a
abs a
xa -> a -> a
forall a. Floating a => a -> a -> a
**a
1.5)
derivative Transformation
Exp     = a -> a
forall a. Floating a => a -> a
exp
derivative Transformation
Log     = a -> a
forall a. Fractional a => a -> a
recip
derivative Transformation
Log1p   = a -> a
forall a. Fractional a => a -> a
recip (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a
forall a. Num a => a -> a -> a
+a
1)

-- | Evaluate the second derivative of the transformation function. 
sndDerivative :: Floating a => Transformation -> a -> a
sndDerivative :: Transformation -> a -> a
sndDerivative Transformation
Id      = a -> a -> a
forall a b. a -> b -> a
const a
0
sndDerivative Transformation
Sin     = a -> a
forall a. Num a => a -> a
negate(a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.a -> a
forall a. Floating a => a -> a
sin
sndDerivative Transformation
Cos     = a -> a
forall a. Num a => a -> a
negate(a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.a -> a
forall a. Floating a => a -> a
cos
sndDerivative Transformation
Tan     = \a
x -> a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
tan a
x a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
cos a
x a -> a -> a
forall a. Floating a => a -> a -> a
** (-a
2.0)
sndDerivative Transformation
Tanh    = \a
x -> (-a
2) a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
tanh a
x 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
tanh a
x a -> a -> a
forall a. Floating a => a -> a -> a
** a
2.0)
sndDerivative Transformation
Sqrt    = a -> a
forall a. Num a => a -> a
negate (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> a
forall a. Fractional a => a -> a
recip (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a
forall a. Num a => a -> a -> a
*a
4) (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> a
forall a. Floating a => a -> a
sqrt (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a
forall a. Floating a => a -> a -> a
**a
3.0)
sndDerivative Transformation
SqrtAbs = \a
x -> -(a
xa -> a -> a
forall a. Floating a => a -> a -> a
**a
2.0 a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
4 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Num a => a -> a
abs a
x a -> a -> a
forall a. Floating a => a -> a -> a
** a
3.5)) -- assuming dirac(x) = 0
sndDerivative Transformation
Exp     = a -> a
forall a. Floating a => a -> a
exp
sndDerivative Transformation
Log     = \a
x -> -(a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
xa -> a -> a
forall a. Floating a => a -> a -> a
**a
2)
sndDerivative Transformation
Log1p   = a -> a
forall a. Num a => a -> a
negate (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> a
forall a. Fractional a => a -> a
recip (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a
forall a. Floating a => a -> a -> a
**a
2.0) (a -> a) -> (a -> a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a
forall a. Num a => a -> a -> a
+a
1)

-- | Evaluates the interaction 'ks' in the dataset 'xss'
--
-- It folds the map of exponents with the key, retrieve the corresponding column,
-- calculates x_i^k_i and multiplies to the accumulator.
-- It return a column vector. 
polynomial :: Dataset Double -> Interaction -> Column Double
polynomial :: Dataset Double -> Interaction -> Column Double
polynomial Dataset Double
xss Interaction
ks
  | Dataset Double -> Int
forall a. Vector a -> Int
V.length Dataset Double
xss Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = [Double] -> Column Double
forall a. Storable a => [a] -> Vector a
LA.fromList []
  | Bool
otherwise         = (Int -> Int -> Column Double -> Column Double)
-> Column Double -> Interaction -> Column Double
forall a b. (Int -> a -> b -> b) -> b -> IntMap a -> b
M.foldrWithKey Int -> Int -> Column Double -> Column Double
forall b. Integral b => Int -> b -> Column Double -> Column Double
monoProduct Column Double
1 Interaction
ks

  where
    monoProduct :: Int -> b -> Column Double -> Column Double
monoProduct Int
ix b
k Column Double
ps = Column Double
ps Column Double -> Column Double -> Column Double
forall a. Num a => a -> a -> a
* (Dataset Double
xss Dataset Double -> Int -> Column Double
forall a. Vector a -> Int -> a
V.! Int
ix) Column Double -> b -> Column Double
forall a b. (Fractional a, Integral b) => a -> b -> a
^^ b
k

-- | Evaluates the interval of the interaction w.r.t. the domains of the variables.
--
-- In this case it is faster to fold through the list of domains and checking whether we 
-- have a nonzero strength.
polynomialInterval :: [Interval Double] -> Interaction -> Interval Double
polynomialInterval :: [Interval Double] -> Interaction -> Interval Double
polynomialInterval [Interval Double]
domains Interaction
ks = ((Int, Interval Double) -> Interval Double -> Interval Double)
-> Interval Double -> [(Int, Interval Double)] -> Interval Double
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (Int, Interval Double) -> Interval Double -> Interval Double
forall p. Floating p => (Int, p) -> p -> p
monoProduct (Double -> Interval Double
forall a. a -> Interval a
singleton Double
1) ([(Int, Interval Double)] -> Interval Double)
-> [(Int, Interval Double)] -> Interval Double
forall a b. (a -> b) -> a -> b
$ [Int] -> [Interval Double] -> [(Int, Interval Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0..] [Interval Double]
domains
  where
    monoProduct :: (Int, p) -> p -> p
monoProduct (Int
ix, p
x) p
img
      | Int
ix Int -> Interaction -> Bool
forall a. Int -> IntMap a -> Bool
`M.member` Interaction
ks = p
img p -> p -> p
forall a. Num a => a -> a -> a
* p
x p -> p -> p
forall a. Floating a => a -> a -> a
** Int -> p
forall b. Num b => Int -> b
get Int
ix
      | Bool
otherwise        = p
img
    get :: Int -> b
get Int
ix = Int -> b
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Interaction
ks Interaction -> Int -> Int
forall a. IntMap a -> Int -> a
M.! Int
ix)

-- | Evaluates a term by applying the transformation function to
-- the result of the interaction.
evalTerm :: Dataset Double -> Term -> Column Double
evalTerm :: Dataset Double -> Term -> Column Double
evalTerm Dataset Double
xss (Term Transformation
t Interaction
ks) = Transformation -> Column Double -> Column Double
forall a. Floating a => Transformation -> a -> a
transform Transformation
t (Dataset Double -> Interaction -> Column Double
polynomial Dataset Double
xss Interaction
ks)

-- | Evaluates the term on the interval domains. 
evalTermInterval :: [Interval Double] -> Term -> Interval Double
evalTermInterval :: [Interval Double] -> Term -> Interval Double
evalTermInterval [Interval Double]
domains (Term Transformation
t Interaction
ks)
  | Interval Double
inter Interval Double -> Interval Double -> Bool
forall a. Eq a => a -> a -> Bool
== Interval Double
forall a. Interval a
empty = Interval Double
forall a. Interval a
empty
  | Bool
otherwise      = (Interval Double -> Interval Double)
-> Interval Double -> Interval Double
forall a.
RealFloat a =>
(Interval a -> Interval a) -> Interval a -> Interval a
protected (Transformation -> Interval Double -> Interval Double
forall a. Floating a => Transformation -> a -> a
transform Transformation
t) Interval Double
inter
  where inter :: Interval Double
inter = [Interval Double] -> Interaction -> Interval Double
polynomialInterval [Interval Double]
domains Interaction
ks

-- | The partial derivative of a term w.r.t. the variable ix is:
-- 
-- If ix exists in the interaction, the derivative is:
--      derivative t (polynomial xss ks) * (polynomial xss ks') 
--      ks' is the same as ks but with the strength of ix decremented by one
-- Otherwise it is equal to 0
--
-- Given the expression: w1 t1(p1(x)) + w2 t2(p2(x))
-- The derivative is: w1 t1'(p1(x))p1'(x) + w2 t2'(p2(x))p2'(x)
evalTermDiff :: Int -> Dataset Double -> Term -> Column Double
evalTermDiff :: Int -> Dataset Double -> Term -> Column Double
evalTermDiff Int
ix Dataset Double
xss (Term Transformation
t Interaction
ks)
  | Int -> Interaction -> Bool
forall a. Int -> IntMap a -> Bool
M.member Int
ix Interaction
ks = Column Double
ki Column Double -> Column Double -> Column Double
forall a. Num a => a -> a -> a
* Column Double
it Column Double -> Column Double -> Column Double
forall a. Num a => a -> a -> a
* Column Double
p'
  | Bool
otherwise      = [Double] -> Column Double
forall a. Storable a => [a] -> Vector a
LA.fromList ([Double] -> Column Double) -> [Double] -> Column Double
forall a b. (a -> b) -> a -> b
$ Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
n Double
0
  where
    ki :: Column Double
ki = Int -> Column Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Column Double) -> Int -> Column Double
forall a b. (a -> b) -> a -> b
$ Interaction
ks Interaction -> Int -> Int
forall a. IntMap a -> Int -> a
M.! Int
ix
    p :: Column Double
p  = Dataset Double -> Interaction -> Column Double
polynomial Dataset Double
xss Interaction
ks
    p' :: Column Double
p' = Dataset Double -> Interaction -> Column Double
polynomial Dataset Double
xss ((Int -> Maybe Int) -> Int -> Interaction -> Interaction
forall a. (a -> Maybe a) -> Int -> IntMap a -> IntMap a
M.update Int -> Maybe Int
forall a. (Eq a, Num a) => a -> Maybe a
dec Int
ix Interaction
ks)
    t' :: Double -> Double
t' = Transformation -> Double -> Double
forall a. Floating a => Transformation -> a -> a
derivative Transformation
t
    it :: Column Double
it = (Double -> Double) -> Column Double -> Column Double
forall b (c :: * -> *) e.
(Element b, Container c e) =>
(e -> b) -> c e -> c b
LA.cmap Double -> Double
t' Column Double
p
    n :: IndexOf Vector
n  = Column Double -> IndexOf Vector
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
LA.size (Dataset Double
xss Dataset Double -> Int -> Column Double
forall a. Vector a -> Int -> a
V.! Int
0)

    dec :: a -> Maybe a
dec a
1 = Maybe a
forall a. Maybe a
Nothing
    dec a
k = a -> Maybe a
forall a. a -> Maybe a
Just (a
ka -> a -> a
forall a. Num a => a -> a -> a
-a
1)

-- | The second partial derivative of a term w.r.t. the variables ix and iy is:
-- 
-- given p(x) as the interaction function, t(x) the transformation function,
-- f'(x)|i the first derivative of a function given i and f''(x) the second derivative. 
--
-- If iy exists in p(x) and ix exists in p'(x)|iy, the derivative is:
--
--     kx*ky*t''(x)*p'(x)|ix * p'(x)|iy + kxy*t'(x)p''(x)
--
-- where kx, ky are the original exponents of variables ix and iy,
--     and kxy is ky * kx', with kx' the exponent of ix on p'(x)|iy
--
evalTermSndDiff :: Int -> Int -> Dataset Double -> Term -> Column Double
evalTermSndDiff :: Int -> Int -> Dataset Double -> Term -> Column Double
evalTermSndDiff Int
ix Int
iy Dataset Double
xss (Term Transformation
t Interaction
ks)
  | Int -> Interaction -> Bool
forall a. Int -> IntMap a -> Bool
M.member Int
ix Interaction
ks' Bool -> Bool -> Bool
&& Int -> Interaction -> Bool
forall a. Int -> IntMap a -> Bool
M.member Int
iy Interaction
ks = Column Double
kiColumn Double -> Column Double -> Column Double
forall a. Num a => a -> a -> a
*Column Double
kjColumn Double -> Column Double -> Column Double
forall a. Num a => a -> a -> a
*Column Double
tp''Column Double -> Column Double -> Column Double
forall a. Num a => a -> a -> a
*Column Double
px'Column Double -> Column Double -> Column Double
forall a. Num a => a -> a -> a
*Column Double
py' Column Double -> Column Double -> Column Double
forall a. Num a => a -> a -> a
+ Column Double
kijColumn Double -> Column Double -> Column Double
forall a. Num a => a -> a -> a
*Column Double
tp'Column Double -> Column Double -> Column Double
forall a. Num a => a -> a -> a
*Column Double
pxy'
  | Bool
otherwise                         = [Double] -> Column Double
forall a. Storable a => [a] -> Vector a
LA.fromList ([Double] -> Column Double) -> [Double] -> Column Double
forall a b. (a -> b) -> a -> b
$ Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
n Double
0
  where
    ks' :: Interaction
ks'  = (Int -> Maybe Int) -> Int -> Interaction -> Interaction
forall a. (a -> Maybe a) -> Int -> IntMap a -> IntMap a
M.update Int -> Maybe Int
forall a. (Eq a, Num a) => a -> Maybe a
dec Int
iy Interaction
ks
    ki :: Column Double
ki = Int -> Column Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Column Double) -> Int -> Column Double
forall a b. (a -> b) -> a -> b
$ Interaction
ks Interaction -> Int -> Int
forall a. IntMap a -> Int -> a
M.! Int
ix
    kj :: Column Double
kj =  Int -> Column Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Column Double) -> Int -> Column Double
forall a b. (a -> b) -> a -> b
$  Interaction
ks Interaction -> Int -> Int
forall a. IntMap a -> Int -> a
M.! Int
iy
    kij :: Column Double
kij = Column Double
kj Column Double -> Column Double -> Column Double
forall a. Num a => a -> a -> a
* Int -> Column Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Interaction
ks' Interaction -> Int -> Int
forall a. IntMap a -> Int -> a
M.! Int
ix)
    p :: Column Double
p    = Dataset Double -> Interaction -> Column Double
polynomial Dataset Double
xss Interaction
ks
    px' :: Column Double
px'  = Dataset Double -> Interaction -> Column Double
polynomial Dataset Double
xss ((Int -> Maybe Int) -> Int -> Interaction -> Interaction
forall a. (a -> Maybe a) -> Int -> IntMap a -> IntMap a
M.update Int -> Maybe Int
forall a. (Eq a, Num a) => a -> Maybe a
dec Int
ix Interaction
ks)
    py' :: Column Double
py'  = Dataset Double -> Interaction -> Column Double
polynomial Dataset Double
xss ((Int -> Maybe Int) -> Int -> Interaction -> Interaction
forall a. (a -> Maybe a) -> Int -> IntMap a -> IntMap a
M.update Int -> Maybe Int
forall a. (Eq a, Num a) => a -> Maybe a
dec Int
iy Interaction
ks)
    pxy' :: Column Double
pxy' = Dataset Double -> Interaction -> Column Double
polynomial Dataset Double
xss ((Int -> Maybe Int) -> Int -> Interaction -> Interaction
forall a. (a -> Maybe a) -> Int -> IntMap a -> IntMap a
M.update Int -> Maybe Int
forall a. (Eq a, Num a) => a -> Maybe a
dec Int
ix Interaction
ks')
    t' :: Double -> Double
t'   = Transformation -> Double -> Double
forall a. Floating a => Transformation -> a -> a
derivative Transformation
t
    t'' :: Double -> Double
t''  = Transformation -> Double -> Double
forall a. Floating a => Transformation -> a -> a
sndDerivative Transformation
t
    tp' :: Column Double
tp'  = (Double -> Double) -> Column Double -> Column Double
forall b (c :: * -> *) e.
(Element b, Container c e) =>
(e -> b) -> c e -> c b
LA.cmap Double -> Double
t'  Column Double
p
    tp'' :: Column Double
tp'' = (Double -> Double) -> Column Double -> Column Double
forall b (c :: * -> *) e.
(Element b, Container c e) =>
(e -> b) -> c e -> c b
LA.cmap Double -> Double
t'' Column Double
p
    n :: IndexOf Vector
n    = Column Double -> IndexOf Vector
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
LA.size (Dataset Double
xss Dataset Double -> Int -> Column Double
forall a. Vector a -> Int -> a
V.! Int
0)

    dec :: a -> Maybe a
dec a
1 = Maybe a
forall a. Maybe a
Nothing
    dec a
k = a -> Maybe a
forall a. a -> Maybe a
Just (a
ka -> a -> a
forall a. Num a => a -> a -> a
-a
1)

-- | Same as 'evalTermDiff' but optimized for intervals. 
evalTermDiffInterval :: Int -> [Interval Double] -> Term -> Interval Double
evalTermDiffInterval :: Int -> [Interval Double] -> Term -> Interval Double
evalTermDiffInterval Int
ix [Interval Double]
domains (Term Transformation
t Interaction
ks)
  | Int -> Interaction -> Bool
forall a. Int -> IntMap a -> Bool
M.member Int
ix Interaction
ks = Interval Double
ki Interval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
* Interval Double
it Interval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
* Interval Double
p'
  | Bool
otherwise      = Double -> Interval Double
forall a. a -> Interval a
singleton Double
0
  where
    ki :: Interval Double
ki = Double -> Interval Double
forall a. a -> Interval a
singleton (Double -> Interval Double)
-> (Int -> Double) -> Int -> Interval Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Interval Double) -> Int -> Interval Double
forall a b. (a -> b) -> a -> b
$ Interaction
ks Interaction -> Int -> Int
forall a. IntMap a -> Int -> a
M.! Int
ix
    p :: Interval Double
p  = [Interval Double] -> Interaction -> Interval Double
polynomialInterval [Interval Double]
domains Interaction
ks
    p' :: Interval Double
p' = [Interval Double] -> Interaction -> Interval Double
polynomialInterval [Interval Double]
domains ((Int -> Maybe Int) -> Int -> Interaction -> Interaction
forall a. (a -> Maybe a) -> Int -> IntMap a -> IntMap a
M.update Int -> Maybe Int
forall a. (Eq a, Num a) => a -> Maybe a
dec Int
ix Interaction
ks)
    t' :: Interval Double -> Interval Double
t' = (Interval Double -> Interval Double)
-> Interval Double -> Interval Double
forall a.
RealFloat a =>
(Interval a -> Interval a) -> Interval a -> Interval a
protected (Transformation -> Interval Double -> Interval Double
forall a. Floating a => Transformation -> a -> a
derivative Transformation
t)
    it :: Interval Double
it = Interval Double -> Interval Double
t' Interval Double
p

    dec :: a -> Maybe a
dec a
1 = Maybe a
forall a. Maybe a
Nothing
    dec a
k = a -> Maybe a
forall a. a -> Maybe a
Just (a
ka -> a -> a
forall a. Num a => a -> a -> a
-a
1)

-- | Same as 'evalTermSndDiff' but optimized for intervals 
evalTermSndDiffInterval :: Int -> Int -> [Interval Double] -> Term -> Interval Double
evalTermSndDiffInterval :: Int -> Int -> [Interval Double] -> Term -> Interval Double
evalTermSndDiffInterval Int
ix Int
iy [Interval Double]
domains (Term Transformation
t Interaction
ks)
  | Int -> Interaction -> Bool
forall a. Int -> IntMap a -> Bool
M.member Int
ix Interaction
ks' Bool -> Bool -> Bool
&& Int -> Interaction -> Bool
forall a. Int -> IntMap a -> Bool
M.member Int
iy Interaction
ks = Interval Double
kiInterval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
*Interval Double
kjInterval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
*Interval Double
tp''Interval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
*Interval Double
px'Interval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
*Interval Double
py' Interval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
+ Interval Double
kijInterval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
*Interval Double
tp'Interval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
*Interval Double
pxy'
  | Bool
otherwise                         = Double -> Interval Double
forall a. a -> Interval a
singleton Double
0
  where
    ks' :: Interaction
ks' = (Int -> Maybe Int) -> Int -> Interaction -> Interaction
forall a. (a -> Maybe a) -> Int -> IntMap a -> IntMap a
M.update Int -> Maybe Int
forall a. (Eq a, Num a) => a -> Maybe a
dec Int
iy Interaction
ks
    ki :: Interval Double
ki = Double -> Interval Double
forall a. a -> Interval a
singleton (Double -> Interval Double)
-> (Int -> Double) -> Int -> Interval Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Interval Double) -> Int -> Interval Double
forall a b. (a -> b) -> a -> b
$ Interaction
ks Interaction -> Int -> Int
forall a. IntMap a -> Int -> a
M.! Int
ix
    kj :: Interval Double
kj = Double -> Interval Double
forall a. a -> Interval a
singleton (Double -> Interval Double)
-> (Int -> Double) -> Int -> Interval Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Interval Double) -> Int -> Interval Double
forall a b. (a -> b) -> a -> b
$ Interaction
ks Interaction -> Int -> Int
forall a. IntMap a -> Int -> a
M.! Int
iy
    kij :: Interval Double
kij = Interval Double
kj Interval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
* Double -> Interval Double
forall a. a -> Interval a
singleton (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Interaction
ks' Interaction -> Int -> Int
forall a. IntMap a -> Int -> a
M.! Int
ix))
    p :: Interval Double
p    = [Interval Double] -> Interaction -> Interval Double
polynomialInterval [Interval Double]
domains Interaction
ks
    px' :: Interval Double
px'  = [Interval Double] -> Interaction -> Interval Double
polynomialInterval [Interval Double]
domains ((Int -> Maybe Int) -> Int -> Interaction -> Interaction
forall a. (a -> Maybe a) -> Int -> IntMap a -> IntMap a
M.update Int -> Maybe Int
forall a. (Eq a, Num a) => a -> Maybe a
dec Int
ix Interaction
ks)
    py' :: Interval Double
py'  = [Interval Double] -> Interaction -> Interval Double
polynomialInterval [Interval Double]
domains ((Int -> Maybe Int) -> Int -> Interaction -> Interaction
forall a. (a -> Maybe a) -> Int -> IntMap a -> IntMap a
M.update Int -> Maybe Int
forall a. (Eq a, Num a) => a -> Maybe a
dec Int
iy Interaction
ks)
    pxy' :: Interval Double
pxy' = [Interval Double] -> Interaction -> Interval Double
polynomialInterval [Interval Double]
domains ((Int -> Maybe Int) -> Int -> Interaction -> Interaction
forall a. (a -> Maybe a) -> Int -> IntMap a -> IntMap a
M.update Int -> Maybe Int
forall a. (Eq a, Num a) => a -> Maybe a
dec Int
ix (Interaction -> Interaction) -> Interaction -> Interaction
forall a b. (a -> b) -> a -> b
$ (Int -> Maybe Int) -> Int -> Interaction -> Interaction
forall a. (a -> Maybe a) -> Int -> IntMap a -> IntMap a
M.update Int -> Maybe Int
forall a. (Eq a, Num a) => a -> Maybe a
dec Int
iy Interaction
ks)
    t' :: Interval Double -> Interval Double
t'   = (Interval Double -> Interval Double)
-> Interval Double -> Interval Double
forall a.
RealFloat a =>
(Interval a -> Interval a) -> Interval a -> Interval a
protected (Transformation -> Interval Double -> Interval Double
forall a. Floating a => Transformation -> a -> a
derivative Transformation
t)
    t'' :: Interval Double -> Interval Double
t''  = (Interval Double -> Interval Double)
-> Interval Double -> Interval Double
forall a.
RealFloat a =>
(Interval a -> Interval a) -> Interval a -> Interval a
protected (Transformation -> Interval Double -> Interval Double
forall a. Floating a => Transformation -> a -> a
sndDerivative Transformation
t)
    tp' :: Interval Double
tp'  = Interval Double -> Interval Double
t'  Interval Double
p
    tp'' :: Interval Double
tp'' = Interval Double -> Interval Double
t'' Interval Double
p

    dec :: a -> Maybe a
dec a
1 = Maybe a
forall a. Maybe a
Nothing
    dec a
k = a -> Maybe a
forall a. a -> Maybe a
Just (a
ka -> a -> a
forall a. Num a => a -> a -> a
-a
1)

-- | evaluates an expression by evaluating the terms into a list
-- applying the weight and summing the results.
evalGeneric :: (Dataset Double -> Term -> Column Double) -> Dataset Double -> Expr -> [Double] -> Column Double
evalGeneric :: (Dataset Double -> Term -> Column Double)
-> Dataset Double -> Expr -> [Double] -> Column Double
evalGeneric Dataset Double -> Term -> Column Double
f Dataset Double
xss Expr
terms [Double]
ws = Column Double
b Column Double -> Column Double -> Column Double
forall a. Num a => a -> a -> a
+ [Column Double] -> Column Double
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Column Double]
weightedTerms
  where
    weightedTerms :: [Column Double]
weightedTerms = (Double -> Column Double -> Column Double)
-> [Double] -> [Column Double] -> [Column Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Column Double -> Column Double
forall b (c :: * -> *). (Container c b, Num b) => b -> c b -> c b
multWeight [Double]
ws ((Term -> Column Double) -> Expr -> [Column Double]
forall a b. (a -> b) -> [a] -> [b]
map (Dataset Double -> Term -> Column Double
f Dataset Double
xss') Expr
terms)
    multWeight :: b -> c b -> c b
multWeight b
w  = (b -> b) -> c b -> c b
forall b (c :: * -> *) e.
(Element b, Container c e) =>
(e -> b) -> c e -> c b
LA.cmap (b
wb -> b -> b
forall a. Num a => a -> a -> a
*)
    b :: Column Double
b             = Dataset Double -> Column 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

-- | Evaluating an expression is simply applying 'evalTerm' 
evalExpr :: Dataset Double -> Expr -> [Double] -> Column Double
evalExpr :: Dataset Double -> Expr -> [Double] -> Column Double
evalExpr = (Dataset Double -> Term -> Column Double)
-> Dataset Double -> Expr -> [Double] -> Column Double
evalGeneric Dataset Double -> Term -> Column Double
evalTerm

-- | Evaluating the derivative is applying 'evalTermDiff' in the expression 
evalDiff :: Int -> Dataset Double -> Expr -> [Double] -> Column Double
evalDiff :: Int -> Dataset Double -> Expr -> [Double] -> Column Double
evalDiff Int
ix = (Dataset Double -> Term -> Column Double)
-> Dataset Double -> Expr -> [Double] -> Column Double
evalGeneric (Int -> Dataset Double -> Term -> Column Double
evalTermDiff Int
ix)

-- | Evaluating the second derivative is applying 'evalTermSndDiff' in the expression
evalSndDiff :: Int -> Int -> Dataset Double -> Expr -> [Double] -> Column Double
evalSndDiff :: Int -> Int -> Dataset Double -> Expr -> [Double] -> Column Double
evalSndDiff Int
ix Int
iy = (Dataset Double -> Term -> Column Double)
-> Dataset Double -> Expr -> [Double] -> Column Double
evalGeneric (Int -> Int -> Dataset Double -> Term -> Column Double
evalTermSndDiff Int
ix Int
iy)

-- | Returns the estimate of the image of the funcion with Interval Arithmetic
--
-- this requires a different implementation from `evalGeneric`
evalImageGeneric :: ([Interval Double] -> Term -> Interval Double) -> [Interval Double]-> Expr -> [Double] -> Interval Double
evalImageGeneric :: ([Interval Double] -> Term -> Interval Double)
-> [Interval Double] -> Expr -> [Double] -> Interval Double
evalImageGeneric [Interval Double] -> Term -> Interval Double
f [Interval Double]
domains Expr
terms [Double]
ws = [Interval Double] -> Interval Double
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Interval Double]
weightedTerms
  where
    weightedTerms :: [Interval Double]
weightedTerms = (Interval Double -> Interval Double -> Interval Double)
-> [Interval Double] -> [Interval Double] -> [Interval Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Interval Double -> Interval Double -> Interval Double
forall a. Num a => a -> a -> a
(*) ((Double -> Interval Double) -> [Double] -> [Interval Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Interval Double
forall a. a -> Interval a
singleton [Double]
ws) ((Term -> Interval Double) -> Expr -> [Interval Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Interval Double] -> Term -> Interval Double
f [Interval Double]
domains) Expr
terms)

evalImage :: [Interval Double] -> Expr -> [Double] -> Interval Double
evalImage :: [Interval Double] -> Expr -> [Double] -> Interval Double
evalImage = ([Interval Double] -> Term -> Interval Double)
-> [Interval Double] -> Expr -> [Double] -> Interval Double
evalImageGeneric [Interval Double] -> Term -> Interval Double
evalTermInterval

evalDiffImage :: Int -> [Interval Double] -> Expr -> [Double] -> Interval Double
evalDiffImage :: Int -> [Interval Double] -> Expr -> [Double] -> Interval Double
evalDiffImage Int
ix = ([Interval Double] -> Term -> Interval Double)
-> [Interval Double] -> Expr -> [Double] -> Interval Double
evalImageGeneric (Int -> [Interval Double] -> Term -> Interval Double
evalTermDiffInterval Int
ix)

evalSndDiffImage :: Int -> Int -> [Interval Double] -> Expr -> [Double] -> Interval Double
evalSndDiffImage :: Int
-> Int -> [Interval Double] -> Expr -> [Double] -> Interval Double
evalSndDiffImage Int
ix Int
iy = ([Interval Double] -> Term -> Interval Double)
-> [Interval Double] -> Expr -> [Double] -> Interval Double
evalImageGeneric (Int -> Int -> [Interval Double] -> Term -> Interval Double
evalTermSndDiffInterval Int
ix Int
iy)

-- | A value is invalid if it's wether NaN or Infinite
isInvalid :: Double -> Bool
isInvalid :: Double -> Bool
isInvalid Double
x = Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x Bool -> Bool -> Bool
|| Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
1e150

-- | a set of points is valid if none of its values are invalid and
-- the maximum abosolute value is below 1e150 (to avoid overflow)
isValid :: [Double] -> Bool
isValid :: [Double] -> Bool
isValid = (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
x -> Bool -> Bool
not (Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x) Bool -> Bool -> Bool
&& Bool -> Bool
not (Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x) Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1e150)
-- not (any isInvalid xs)

-- | evaluate an expression to a set of samples 
--
-- (1 LA.|||) adds a bias dimension
exprToMatrix :: Dataset Double -> Expr -> LA.Matrix Double
exprToMatrix :: Dataset Double -> Expr -> Matrix Double
exprToMatrix Dataset Double
xss = [Column Double] -> Matrix Double
forall t. Element t => [Vector t] -> Matrix t
LA.fromColumns ([Column Double] -> Matrix Double)
-> (Expr -> [Column Double]) -> Expr -> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Dataset Double -> Column Double
forall a. Vector a -> a
V.head Dataset Double
xss Column Double -> [Column Double] -> [Column Double]
forall a. a -> [a] -> [a]
:) ([Column Double] -> [Column Double])
-> (Expr -> [Column Double]) -> Expr -> [Column Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Term -> Column Double) -> Expr -> [Column Double]
forall a b. (a -> b) -> [a] -> [b]
map (Dataset Double -> Term -> Column Double
evalTerm (Dataset Double -> Dataset Double
forall a. Vector a -> Vector a
V.tail Dataset Double
xss))

-- | Clean the expression by removing the invalid terms
-- TODO: move to IT.Regression 
cleanExpr :: Dataset Double -> Expr -> (Expr, LA.Matrix Double)
cleanExpr :: Dataset Double -> Expr -> (Expr, Matrix Double)
cleanExpr Dataset Double
xss = ([Column Double] -> Matrix Double)
-> (Expr, [Column Double]) -> (Expr, Matrix Double)
forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second ([Column Double] -> Matrix Double
forall t. Element t => [Vector t] -> Matrix t
LA.fromColumns ([Column Double] -> Matrix Double)
-> ([Column Double] -> [Column Double])
-> [Column Double]
-> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Column Double
bColumn Double -> [Column Double] -> [Column Double]
forall a. a -> [a] -> [a]
:)) ((Expr, [Column Double]) -> (Expr, Matrix Double))
-> (Expr -> (Expr, [Column Double]))
-> Expr
-> (Expr, Matrix Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Term -> (Expr, [Column Double]) -> (Expr, [Column Double]))
-> (Expr, [Column Double]) -> Expr -> (Expr, [Column Double])
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr Term -> (Expr, [Column Double]) -> (Expr, [Column Double])
p ([], [])
  where
    xss' :: Dataset Double
xss'           = Dataset Double -> Dataset Double
forall a. Vector a -> Vector a
V.tail Dataset Double
xss
    b :: Column Double
b              = Dataset Double -> Column Double
forall a. Vector a -> a
V.head Dataset Double
xss
    p :: Term -> (Expr, [Column Double]) -> (Expr, [Column Double])
p Term
t (Expr
ts, [Column Double]
cols) = let col :: Column Double
col = Dataset Double -> Term -> Column Double
evalTerm Dataset Double
xss' Term
t
                     in  if (Double -> Bool) -> Column Double -> Bool
forall a. Storable a => (a -> Bool) -> Vector a -> Bool
VV.all (Bool -> Bool
not(Bool -> Bool) -> (Double -> Bool) -> Double -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Double -> Bool
isInvalid) Column Double
col
                            then (Term
tTerm -> Expr -> Expr
forall a. a -> [a] -> [a]
:Expr
ts, Column Double
colColumn Double -> [Column Double] -> [Column Double]
forall a. a -> [a] -> [a]
:[Column Double]
cols)
                            else (Expr
ts, [Column Double]
cols)

-- | Checks if the fitness of a solution is not Inf nor NaN.
notInfNan :: Solution -> Bool
notInfNan :: Solution -> Bool
notInfNan Solution
s = Bool -> Bool
not (Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
f Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
f)
  where f :: Double
f = [Double] -> Double
forall a. [a] -> a
head ([Double] -> Double) -> [Double] -> Double
forall a b. (a -> b) -> a -> b
$ Solution -> [Double]
_fit Solution
s

-- | definition of an interval evaluated to NaN
nanInterval :: RealFloat a => Interval a
nanInterval :: Interval a
nanInterval = a -> Interval a
forall a. a -> Interval a
singleton (a
0a -> a -> a
forall a. Fractional a => a -> a -> a
/a
0)

-- | Check if any bound of the interval is infinity
hasInf :: RealFloat a => Interval a -> Bool
hasInf :: Interval a -> Bool
hasInf Interval a
x = (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite [Interval a -> a
forall a. Interval a -> a
inf Interval a
x, Interval a -> a
forall a. Interval a -> a
sup Interval a
x]

-- | Creates a protected function to avoid throwing exceptions
protected :: RealFloat a => (Interval a -> Interval a) -> Interval a -> Interval a
protected :: (Interval a -> Interval a) -> Interval a -> Interval a
protected Interval a -> Interval a
f Interval a
x
  | Interval a -> Bool
forall a. RealFloat a => Interval a -> Bool
hasInf Interval a
x      = Interval a
forall a. RealFloat a => Interval a
nanInterval
  | Interval a -> a
forall a. Interval a -> a
sup Interval a
y a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< Interval a -> a
forall a. Interval a -> a
inf Interval a
y = Interval a -> a
forall a. Interval a -> a
sup Interval a
y a -> a -> Interval a
forall a. Ord a => a -> a -> Interval a
... Interval a -> a
forall a. Interval a -> a
inf Interval a
y
  | Bool
otherwise     = Interval a
y
  where
    y :: Interval a
y = Interval a -> Interval a
f Interval a
x