```{- | 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
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))

```