{-# OPTIONS_GHC -fglasgow-exts #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.GSL.Fitting.Linear
-- Copyright   :  (c) Alexander Vivian Hugh McPhail 2010
-- License     :  GPL-style
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  uses ffi
--
-- GSL linear regression functions
--
-----------------------------------------------------------------------------

module Numeric.GSL.Fitting.Linear (
                                   linear, linear_w, linear_est
                             ) where

-----------------------------------------------------------------------------

import Data.Packed.Vector
--import Data.Packed.Matrix
import Data.Packed.Development

--import Numeric.LinearAlgebra.Linear

--import Control.Monad(when)

import Foreign
--import Foreign.ForeignPtr
--import Foreign.Marshal.Alloc(alloca)
import Foreign.C.Types(CInt)
--import Foreign.C.String(newCString,peekCString)

--import GHC.ForeignPtr           (mallocPlainForeignPtrBytes)

--import GHC.Base
--import GHC.IOBase

--import Prelude hiding(reverse)

-----------------------------------------------------------------------------

-- | fits the model Y = C X
linear :: Vector Double    -- ^ x data
       -> Vector Double    -- ^ y data
       -> (Double,Double,Double,Double,Double,Double) -- ^ (c_0,c_1,cov_00,cov_01,cov_11,chi_sq)
linear x y = unsafePerformIO $ do
           alloca $ \c0 ->
               alloca $ \c1 ->
                   alloca $ \chi_sq -> 
                       alloca $ \cov00 -> 
                           alloca $ \cov01 -> 
                               alloca $ \cov11 -> do
                                                  app2 (fitting_linear c0 c1 chi_sq cov00 cov01 cov11) vec x vec y "linear"
                                                  c0' <- peek c0
                                                  c1' <- peek c1
                                                  cov00' <- peek cov00
                                                  cov01' <- peek cov01
                                                  cov11' <- peek cov11
                                                  chi_sq' <- peek chi_sq
                                                  return (c0',c1',cov00',cov01',cov11',chi_sq')


-----------------------------------------------------------------------------

foreign import ccall "fitting-aux.h linear" fitting_linear :: Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt

-----------------------------------------------------------------------------

-- | fits the model Y = C X, with x data weighted
linear_w :: Vector Double    -- ^ x data
         -> Vector Double    -- ^ x weights
         -> Vector Double    -- ^ y data
         -> (Double,Double,Double,Double,Double,Double) -- ^ (c_0,c_1,cov_00,cov_01,cov_11,chi_sq)
linear_w x w y = unsafePerformIO $ do
             alloca $ \c0 ->
                 alloca $ \c1 ->
                     alloca $ \chi_sq ->
                       alloca $ \cov00 -> 
                           alloca $ \cov01 -> 
                               alloca $ \cov11 -> do
                                                  app3 (fitting_linear_w c0 c1 chi_sq cov00 cov01 cov11) vec x vec w vec y "linear_w"
                                                  c0' <- peek c0
                                                  c1' <- peek c1
                                                  cov00' <- peek cov00
                                                  cov01' <- peek cov01
                                                  cov11' <- peek cov11
                                                  chi_sq' <- peek chi_sq
                                                  return (c0',c1',cov00',cov01',cov11',chi_sq')

-----------------------------------------------------------------------------

foreign import ccall "fitting-aux.h linear_weighted" fitting_linear_w :: Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt

-----------------------------------------------------------------------------

-- | computes the fitted function and standard deviation at the input point
linear_est :: Double          -- ^ x data point
           -> Double          -- ^ c0
           -> Double          -- ^ c1
           -> Double          -- ^ cov00
           -> Double          -- ^ cov01
           -> Double          -- ^ cov11
           -> (Double,Double) -- ^ (y,error)
linear_est x c0 c1 cov00 cov01 cov11 = unsafePerformIO $ do
                                       alloca $ \y ->
                                           alloca $ \e -> do
                                               check "linear_est" $ fitting_linear_est x c0 c1 cov00 cov01 cov11 y e
                                               y' <- peek y
                                               e' <- peek e
                                               return (y',e')

-----------------------------------------------------------------------------

foreign import ccall "fitting-aux.h linear_estimate" fitting_linear_est :: Double -> Double -> Double -> Double -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt

-----------------------------------------------------------------------------