{- | HasGP Gaussian Process Library. This module contains assorted functions 
     that support GP calculations but are more general-purpose than 
     GP-specific.

     Copyright (C) 2011 Sean Holden. sbh11\@cl.cam.ac.uk.
-}
{- This file is part of HasGP.

   HasGP is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   HasGP is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with HasGP.  If not, see <http://www.gnu.org/licenses/>.
-}
module HasGP.Support.Functions where

import HasGP.Types.MainTypes
import Numeric.LinearAlgebra
import Numeric.GSL.Special.Erf

square :: Double -> Double
square x = (x * x)

trace :: DMatrix -> Double
trace = sum . toList . takeDiag

-- | Standard delta function - 0/1 valued.
delta :: (Eq a) => a -> a -> Double
delta a b 
    | (a==b) = 1.0
    | otherwise = 0.0

-- | Standard delta function - boolean valued.
deltaBool :: (Eq a) => a -> a -> Bool
deltaBool a b = (a==b)

-- | General sigmoid function with variable slope.
generalSigmoid :: Double -> Double -> Double
generalSigmoid theta x = 1 / (1 + (exp (-(theta * x))))

-- | Standard sigmoid function.
sigmoid :: Double -> Double
sigmoid = generalSigmoid 1

-- | Integral of Gaussian density of mean 0 and variance 1 
--   from -infinity to x
phiIntegral :: Double -> Double
phiIntegral x = 1 - (erf_Q x)

-- | Value of Gaussian density function for mean 0 and 
--   variance 1.
n :: Double -> Double
n x = erf_Z x

-- | DANGER! You can't compute the ratio (n x) / (phiIntegral x) directly, 
--   as although it has sensible values for negative x the denominator gets 
--   small so fast that you quickly get Infinity turning up. GSL has the 
--   inverse Mill's function/hazard function for the Gaussian distribution, 
--   and the ratio is equal to hazard(-x).
nOverPhi :: Double -> Double
nOverPhi x = hazard(-x)

-- | DANGER! See nOverPhi - you have to compute this carefully as 
--   well.
logPhi :: Double -> Double
logPhi x = log $ (n x) / (hazard (-x))