{- |
Module      :  Numeric.GSL.Differentiation
Copyright   :  (c) Alberto Ruiz 2006
License     :  GPL

Maintainer  :  Alberto Ruiz
Stability   :  provisional

Numerical differentiation.

<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Differentiation.html#Numerical-Differentiation>

From the GSL manual: \"The functions described in this chapter compute numerical derivatives by finite differencing. An adaptive algorithm is used to find the best choice of finite difference and to estimate the error in the derivative.\"
-}


module Numeric.GSL.Differentiation (
    derivCentral,
    derivForward,
    derivBackward
) where

import Foreign.C.Types
import Foreign.Marshal.Alloc(malloc, free)
import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
import Foreign.Storable(peek)
import Numeric.GSL.Internal
import System.IO.Unsafe(unsafePerformIO)

derivGen ::
    CInt                   -- ^ type: 0 central, 1 forward, 2 backward
    -> Double             -- ^ initial step size
    -> (Double -> Double) -- ^ function
    -> Double             -- ^ point where the derivative is taken
    -> (Double, Double)   -- ^ result and error
derivGen :: CInt -> Double -> (Double -> Double) -> Double -> (Double, Double)
derivGen CInt
c Double
h Double -> Double
f Double
x = IO (Double, Double) -> (Double, Double)
forall a. IO a -> a
unsafePerformIO (IO (Double, Double) -> (Double, Double))
-> IO (Double, Double) -> (Double, Double)
forall a b. (a -> b) -> a -> b
$ do
    Ptr Double
r <- IO (Ptr Double)
forall a. Storable a => IO (Ptr a)
malloc
    Ptr Double
e <- IO (Ptr Double)
forall a. Storable a => IO (Ptr a)
malloc
    FunPtr (Double -> Ptr () -> Double)
fp <- (Double -> Ptr () -> Double)
-> IO (FunPtr (Double -> Ptr () -> Double))
mkfun (\Double
y Ptr ()
_ -> Double -> Double
f Double
y) 
    CInt
-> FunPtr (Double -> Ptr () -> Double)
-> Double
-> Double
-> Ptr Double
-> Ptr Double
-> IO CInt
c_deriv CInt
c FunPtr (Double -> Ptr () -> Double)
fp Double
x Double
h Ptr Double
r Ptr Double
e IO CInt -> (IO CInt -> IO ()) -> IO ()
forall x y. x -> (x -> y) -> y
// String -> IO CInt -> IO ()
check String
"deriv"
    Double
vr <- Ptr Double -> IO Double
forall a. Storable a => Ptr a -> IO a
peek Ptr Double
r
    Double
ve <- Ptr Double -> IO Double
forall a. Storable a => Ptr a -> IO a
peek Ptr Double
e
    let result :: (Double, Double)
result = (Double
vr,Double
ve)
    Ptr Double -> IO ()
forall a. Ptr a -> IO ()
free Ptr Double
r
    Ptr Double -> IO ()
forall a. Ptr a -> IO ()
free Ptr Double
e
    FunPtr (Double -> Ptr () -> Double) -> IO ()
forall a. FunPtr a -> IO ()
freeHaskellFunPtr FunPtr (Double -> Ptr () -> Double)
fp
    (Double, Double) -> IO (Double, Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Double, Double)
result

foreign import ccall safe "gsl-aux.h deriv" 
 c_deriv :: CInt -> FunPtr (Double -> Ptr () -> Double) -> Double -> Double 
                    -> Ptr Double -> Ptr Double -> IO CInt


{- | Adaptive central difference algorithm, /gsl_deriv_central/. For example:

>>> let deriv = derivCentral 0.01
>>> deriv sin (pi/4)
(0.7071067812000676,1.0600063101654055e-10)
>>> cos (pi/4)
0.7071067811865476

-}
derivCentral :: Double                  -- ^ initial step size
                -> (Double -> Double)   -- ^ function 
                -> Double               -- ^ point where the derivative is taken
                -> (Double, Double)     -- ^ result and absolute error
derivCentral :: Double -> (Double -> Double) -> Double -> (Double, Double)
derivCentral = CInt -> Double -> (Double -> Double) -> Double -> (Double, Double)
derivGen CInt
0

-- | Adaptive forward difference algorithm, /gsl_deriv_forward/. The function is evaluated only at points greater than x, and never at x itself. The derivative is returned in result and an estimate of its absolute error is returned in abserr. This function should be used if f(x) has a discontinuity at x, or is undefined for values less than x. A backward derivative can be obtained using a negative step.
derivForward :: Double                  -- ^ initial step size
                -> (Double -> Double)   -- ^ function 
                -> Double               -- ^ point where the derivative is taken
                -> (Double, Double)     -- ^ result and absolute error
derivForward :: Double -> (Double -> Double) -> Double -> (Double, Double)
derivForward = CInt -> Double -> (Double -> Double) -> Double -> (Double, Double)
derivGen CInt
1

-- | Adaptive backward difference algorithm, /gsl_deriv_backward/. 
derivBackward ::Double                  -- ^ initial step size
                -> (Double -> Double)   -- ^ function 
                -> Double               -- ^ point where the derivative is taken
                -> (Double, Double)     -- ^ result and absolute error
derivBackward :: Double -> (Double -> Double) -> Double -> (Double, Double)
derivBackward = CInt -> Double -> (Double -> Double) -> Double -> (Double, Double)
derivGen CInt
2

{- | conversion of Haskell functions into function pointers that can be used in the C side
-}
foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double))