{-# LANGUAGE  FlexibleContexts #-}

{-# OPTIONS_GHC -fno-warn-missing-signatures #-}

{- |
Module      :  Numeric.GSL.Root
Copyright   :  (c) Alberto Ruiz 2009
License     :  GPL
Maintainer  :  Alberto Ruiz
Stability   :  provisional

Multidimensional root finding.

<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html>

The example in the GSL manual:

>>> let rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
>>> let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
>>> sol
[1.0,1.0]
>>> disp 3 path
11x5
 1.000  -10.000  -5.000  11.000  -1050.000
 2.000   -3.976  24.827   4.976     90.203
 3.000   -3.976  24.827   4.976     90.203
 4.000   -3.976  24.827   4.976     90.203
 5.000   -1.274  -5.680   2.274    -73.018
 6.000   -1.274  -5.680   2.274    -73.018
 7.000    0.249   0.298   0.751      2.359
 8.000    0.249   0.298   0.751      2.359
 9.000    1.000   0.878  -0.000     -1.218
10.000    1.000   0.989  -0.000     -0.108
11.000    1.000   1.000   0.000      0.000

-}
-----------------------------------------------------------------------------

module Numeric.GSL.Root (
    uniRoot, UniRootMethod(..),
    uniRootJ, UniRootMethodJ(..),
    root, RootMethod(..),
    rootJ, RootMethodJ(..),
) where

import Numeric.LinearAlgebra.HMatrix
import Numeric.GSL.Internal
import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
import Foreign.C.Types
import System.IO.Unsafe(unsafePerformIO)

-------------------------------------------------------------------------
type TVV = TV (TV Res)
type TVM = TV (TM Res)


data UniRootMethod = Bisection
                   | FalsePos
                   | Brent
                   deriving (Int -> UniRootMethod
UniRootMethod -> Int
UniRootMethod -> [UniRootMethod]
UniRootMethod -> UniRootMethod
UniRootMethod -> UniRootMethod -> [UniRootMethod]
UniRootMethod -> UniRootMethod -> UniRootMethod -> [UniRootMethod]
(UniRootMethod -> UniRootMethod)
-> (UniRootMethod -> UniRootMethod)
-> (Int -> UniRootMethod)
-> (UniRootMethod -> Int)
-> (UniRootMethod -> [UniRootMethod])
-> (UniRootMethod -> UniRootMethod -> [UniRootMethod])
-> (UniRootMethod -> UniRootMethod -> [UniRootMethod])
-> (UniRootMethod
    -> UniRootMethod -> UniRootMethod -> [UniRootMethod])
-> Enum UniRootMethod
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
enumFromThenTo :: UniRootMethod -> UniRootMethod -> UniRootMethod -> [UniRootMethod]
$cenumFromThenTo :: UniRootMethod -> UniRootMethod -> UniRootMethod -> [UniRootMethod]
enumFromTo :: UniRootMethod -> UniRootMethod -> [UniRootMethod]
$cenumFromTo :: UniRootMethod -> UniRootMethod -> [UniRootMethod]
enumFromThen :: UniRootMethod -> UniRootMethod -> [UniRootMethod]
$cenumFromThen :: UniRootMethod -> UniRootMethod -> [UniRootMethod]
enumFrom :: UniRootMethod -> [UniRootMethod]
$cenumFrom :: UniRootMethod -> [UniRootMethod]
fromEnum :: UniRootMethod -> Int
$cfromEnum :: UniRootMethod -> Int
toEnum :: Int -> UniRootMethod
$ctoEnum :: Int -> UniRootMethod
pred :: UniRootMethod -> UniRootMethod
$cpred :: UniRootMethod -> UniRootMethod
succ :: UniRootMethod -> UniRootMethod
$csucc :: UniRootMethod -> UniRootMethod
Enum, UniRootMethod -> UniRootMethod -> Bool
(UniRootMethod -> UniRootMethod -> Bool)
-> (UniRootMethod -> UniRootMethod -> Bool) -> Eq UniRootMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: UniRootMethod -> UniRootMethod -> Bool
$c/= :: UniRootMethod -> UniRootMethod -> Bool
== :: UniRootMethod -> UniRootMethod -> Bool
$c== :: UniRootMethod -> UniRootMethod -> Bool
Eq, Int -> UniRootMethod -> ShowS
[UniRootMethod] -> ShowS
UniRootMethod -> String
(Int -> UniRootMethod -> ShowS)
-> (UniRootMethod -> String)
-> ([UniRootMethod] -> ShowS)
-> Show UniRootMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [UniRootMethod] -> ShowS
$cshowList :: [UniRootMethod] -> ShowS
show :: UniRootMethod -> String
$cshow :: UniRootMethod -> String
showsPrec :: Int -> UniRootMethod -> ShowS
$cshowsPrec :: Int -> UniRootMethod -> ShowS
Show, UniRootMethod
UniRootMethod -> UniRootMethod -> Bounded UniRootMethod
forall a. a -> a -> Bounded a
maxBound :: UniRootMethod
$cmaxBound :: UniRootMethod
minBound :: UniRootMethod
$cminBound :: UniRootMethod
Bounded)

uniRoot :: UniRootMethod
        -> Double
        -> Int
        -> (Double -> Double)
        -> Double
        -> Double
        -> (Double, Matrix Double)
uniRoot :: UniRootMethod
-> Double
-> Int
-> (Double -> Double)
-> Double
-> Double
-> (Double, Matrix Double)
uniRoot method :: UniRootMethod
method epsrel :: Double
epsrel maxit :: Int
maxit fun :: Double -> Double
fun xl :: Double
xl xu :: Double
xu = CInt
-> (Double -> Double)
-> Double
-> Double
-> Double
-> Int
-> (Double, Matrix Double)
uniRootGen (Int -> CInt
fi (UniRootMethod -> Int
forall a. Enum a => a -> Int
fromEnum UniRootMethod
method)) Double -> Double
fun Double
xl Double
xu Double
epsrel Int
maxit

uniRootGen :: CInt
-> (Double -> Double)
-> Double
-> Double
-> Double
-> Int
-> (Double, Matrix Double)
uniRootGen m :: CInt
m f :: Double -> Double
f xl :: Double
xl xu :: Double
xu epsrel :: Double
epsrel maxit :: Int
maxit = IO (Double, Matrix Double) -> (Double, Matrix Double)
forall a. IO a -> a
unsafePerformIO (IO (Double, Matrix Double) -> (Double, Matrix Double))
-> IO (Double, Matrix Double) -> (Double, Matrix Double)
forall a b. (a -> b) -> a -> b
$ do
    FunPtr (Double -> Double)
fp <- (Double -> Double) -> IO (FunPtr (Double -> Double))
mkDoublefun Double -> Double
f
    Matrix Double
rawpath <- Int
-> Int
-> (CInt -> CInt -> Ptr Double -> IO CInt)
-> String
-> IO (Matrix Double)
forall a.
Storable a =>
Int
-> Int
-> (CInt -> CInt -> Ptr a -> IO CInt)
-> String
-> IO (Matrix a)
createMIO Int
maxit 4
                         (CInt
-> FunPtr (Double -> Double)
-> Double
-> CInt
-> Double
-> Double
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
c_root CInt
m FunPtr (Double -> Double)
fp Double
epsrel (Int -> CInt
fi Int
maxit) Double
xl Double
xu)
                         "root"
    let it :: Int
it = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Matrix Double
rawpath Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (Int
maxitInt -> Int -> Int
forall a. Num a => a -> a -> a
-1,0))
        path :: Matrix Double
path = Int -> Matrix Double -> Matrix Double
forall t. Element t => Int -> Matrix t -> Matrix t
takeRows Int
it Matrix Double
rawpath
        [sol :: [Double]
sol] = Matrix Double -> [[Double]]
forall t. Element t => Matrix t -> [[t]]
toLists (Matrix Double -> [[Double]]) -> Matrix Double -> [[Double]]
forall a b. (a -> b) -> a -> b
$ Int -> Matrix Double -> Matrix Double
forall t. Element t => Int -> Matrix t -> Matrix t
dropRows (Int
itInt -> Int -> Int
forall a. Num a => a -> a -> a
-1) Matrix Double
path
    FunPtr (Double -> Double) -> IO ()
forall a. FunPtr a -> IO ()
freeHaskellFunPtr FunPtr (Double -> Double)
fp
    (Double, Matrix Double) -> IO (Double, Matrix Double)
forall (m :: * -> *) a. Monad m => a -> m a
return ([Double]
sol [Double] -> Int -> Double
forall a. [a] -> Int -> a
!! 1, Matrix Double
path)

foreign import ccall safe "root"
    c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM Res

-------------------------------------------------------------------------
data UniRootMethodJ = UNewton
                    | Secant
                    | Steffenson
                    deriving (Int -> UniRootMethodJ
UniRootMethodJ -> Int
UniRootMethodJ -> [UniRootMethodJ]
UniRootMethodJ -> UniRootMethodJ
UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
UniRootMethodJ
-> UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
(UniRootMethodJ -> UniRootMethodJ)
-> (UniRootMethodJ -> UniRootMethodJ)
-> (Int -> UniRootMethodJ)
-> (UniRootMethodJ -> Int)
-> (UniRootMethodJ -> [UniRootMethodJ])
-> (UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ])
-> (UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ])
-> (UniRootMethodJ
    -> UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ])
-> Enum UniRootMethodJ
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
enumFromThenTo :: UniRootMethodJ
-> UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
$cenumFromThenTo :: UniRootMethodJ
-> UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
enumFromTo :: UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
$cenumFromTo :: UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
enumFromThen :: UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
$cenumFromThen :: UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
enumFrom :: UniRootMethodJ -> [UniRootMethodJ]
$cenumFrom :: UniRootMethodJ -> [UniRootMethodJ]
fromEnum :: UniRootMethodJ -> Int
$cfromEnum :: UniRootMethodJ -> Int
toEnum :: Int -> UniRootMethodJ
$ctoEnum :: Int -> UniRootMethodJ
pred :: UniRootMethodJ -> UniRootMethodJ
$cpred :: UniRootMethodJ -> UniRootMethodJ
succ :: UniRootMethodJ -> UniRootMethodJ
$csucc :: UniRootMethodJ -> UniRootMethodJ
Enum, UniRootMethodJ -> UniRootMethodJ -> Bool
(UniRootMethodJ -> UniRootMethodJ -> Bool)
-> (UniRootMethodJ -> UniRootMethodJ -> Bool) -> Eq UniRootMethodJ
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: UniRootMethodJ -> UniRootMethodJ -> Bool
$c/= :: UniRootMethodJ -> UniRootMethodJ -> Bool
== :: UniRootMethodJ -> UniRootMethodJ -> Bool
$c== :: UniRootMethodJ -> UniRootMethodJ -> Bool
Eq, Int -> UniRootMethodJ -> ShowS
[UniRootMethodJ] -> ShowS
UniRootMethodJ -> String
(Int -> UniRootMethodJ -> ShowS)
-> (UniRootMethodJ -> String)
-> ([UniRootMethodJ] -> ShowS)
-> Show UniRootMethodJ
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [UniRootMethodJ] -> ShowS
$cshowList :: [UniRootMethodJ] -> ShowS
show :: UniRootMethodJ -> String
$cshow :: UniRootMethodJ -> String
showsPrec :: Int -> UniRootMethodJ -> ShowS
$cshowsPrec :: Int -> UniRootMethodJ -> ShowS
Show, UniRootMethodJ
UniRootMethodJ -> UniRootMethodJ -> Bounded UniRootMethodJ
forall a. a -> a -> Bounded a
maxBound :: UniRootMethodJ
$cmaxBound :: UniRootMethodJ
minBound :: UniRootMethodJ
$cminBound :: UniRootMethodJ
Bounded)

uniRootJ :: UniRootMethodJ
        -> Double
        -> Int
        -> (Double -> Double)
        -> (Double -> Double)
        -> Double
        -> (Double, Matrix Double)
uniRootJ :: UniRootMethodJ
-> Double
-> Int
-> (Double -> Double)
-> (Double -> Double)
-> Double
-> (Double, Matrix Double)
uniRootJ method :: UniRootMethodJ
method epsrel :: Double
epsrel maxit :: Int
maxit fun :: Double -> Double
fun dfun :: Double -> Double
dfun x :: Double
x = CInt
-> (Double -> Double)
-> (Double -> Double)
-> Double
-> Double
-> Int
-> (Double, Matrix Double)
uniRootJGen (Int -> CInt
fi (UniRootMethodJ -> Int
forall a. Enum a => a -> Int
fromEnum UniRootMethodJ
method)) Double -> Double
fun
    Double -> Double
dfun Double
x Double
epsrel Int
maxit

uniRootJGen :: CInt
-> (Double -> Double)
-> (Double -> Double)
-> Double
-> Double
-> Int
-> (Double, Matrix Double)
uniRootJGen m :: CInt
m f :: Double -> Double
f df :: Double -> Double
df x :: Double
x epsrel :: Double
epsrel maxit :: Int
maxit = IO (Double, Matrix Double) -> (Double, Matrix Double)
forall a. IO a -> a
unsafePerformIO (IO (Double, Matrix Double) -> (Double, Matrix Double))
-> IO (Double, Matrix Double) -> (Double, Matrix Double)
forall a b. (a -> b) -> a -> b
$ do
    FunPtr (Double -> Double)
fp <- (Double -> Double) -> IO (FunPtr (Double -> Double))
mkDoublefun Double -> Double
f
    FunPtr (Double -> Double)
dfp <- (Double -> Double) -> IO (FunPtr (Double -> Double))
mkDoublefun Double -> Double
df
    Matrix Double
rawpath <- Int
-> Int
-> (CInt -> CInt -> Ptr Double -> IO CInt)
-> String
-> IO (Matrix Double)
forall a.
Storable a =>
Int
-> Int
-> (CInt -> CInt -> Ptr a -> IO CInt)
-> String
-> IO (Matrix a)
createMIO Int
maxit 2
                         (CInt
-> FunPtr (Double -> Double)
-> FunPtr (Double -> Double)
-> Double
-> CInt
-> Double
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
c_rootj CInt
m FunPtr (Double -> Double)
fp FunPtr (Double -> Double)
dfp Double
epsrel (Int -> CInt
fi Int
maxit) Double
x)
                         "rootj"
    let it :: Int
it = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Matrix Double
rawpath Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (Int
maxitInt -> Int -> Int
forall a. Num a => a -> a -> a
-1,0))
        path :: Matrix Double
path = Int -> Matrix Double -> Matrix Double
forall t. Element t => Int -> Matrix t -> Matrix t
takeRows Int
it Matrix Double
rawpath
        [sol :: [Double]
sol] = Matrix Double -> [[Double]]
forall t. Element t => Matrix t -> [[t]]
toLists (Matrix Double -> [[Double]]) -> Matrix Double -> [[Double]]
forall a b. (a -> b) -> a -> b
$ Int -> Matrix Double -> Matrix Double
forall t. Element t => Int -> Matrix t -> Matrix t
dropRows (Int
itInt -> Int -> Int
forall a. Num a => a -> a -> a
-1) Matrix Double
path
    FunPtr (Double -> Double) -> IO ()
forall a. FunPtr a -> IO ()
freeHaskellFunPtr FunPtr (Double -> Double)
fp
    (Double, Matrix Double) -> IO (Double, Matrix Double)
forall (m :: * -> *) a. Monad m => a -> m a
return ([Double]
sol [Double] -> Int -> Double
forall a. [a] -> Int -> a
!! 1, Matrix Double
path)

foreign import ccall safe "rootj"
    c_rootj :: CInt -> FunPtr (Double -> Double) -> FunPtr (Double -> Double)
            -> Double -> CInt -> Double -> TM Res

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

data RootMethod = Hybrids
                | Hybrid
                | DNewton
                | Broyden
                deriving (Int -> RootMethod
RootMethod -> Int
RootMethod -> [RootMethod]
RootMethod -> RootMethod
RootMethod -> RootMethod -> [RootMethod]
RootMethod -> RootMethod -> RootMethod -> [RootMethod]
(RootMethod -> RootMethod)
-> (RootMethod -> RootMethod)
-> (Int -> RootMethod)
-> (RootMethod -> Int)
-> (RootMethod -> [RootMethod])
-> (RootMethod -> RootMethod -> [RootMethod])
-> (RootMethod -> RootMethod -> [RootMethod])
-> (RootMethod -> RootMethod -> RootMethod -> [RootMethod])
-> Enum RootMethod
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
enumFromThenTo :: RootMethod -> RootMethod -> RootMethod -> [RootMethod]
$cenumFromThenTo :: RootMethod -> RootMethod -> RootMethod -> [RootMethod]
enumFromTo :: RootMethod -> RootMethod -> [RootMethod]
$cenumFromTo :: RootMethod -> RootMethod -> [RootMethod]
enumFromThen :: RootMethod -> RootMethod -> [RootMethod]
$cenumFromThen :: RootMethod -> RootMethod -> [RootMethod]
enumFrom :: RootMethod -> [RootMethod]
$cenumFrom :: RootMethod -> [RootMethod]
fromEnum :: RootMethod -> Int
$cfromEnum :: RootMethod -> Int
toEnum :: Int -> RootMethod
$ctoEnum :: Int -> RootMethod
pred :: RootMethod -> RootMethod
$cpred :: RootMethod -> RootMethod
succ :: RootMethod -> RootMethod
$csucc :: RootMethod -> RootMethod
Enum,RootMethod -> RootMethod -> Bool
(RootMethod -> RootMethod -> Bool)
-> (RootMethod -> RootMethod -> Bool) -> Eq RootMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: RootMethod -> RootMethod -> Bool
$c/= :: RootMethod -> RootMethod -> Bool
== :: RootMethod -> RootMethod -> Bool
$c== :: RootMethod -> RootMethod -> Bool
Eq,Int -> RootMethod -> ShowS
[RootMethod] -> ShowS
RootMethod -> String
(Int -> RootMethod -> ShowS)
-> (RootMethod -> String)
-> ([RootMethod] -> ShowS)
-> Show RootMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [RootMethod] -> ShowS
$cshowList :: [RootMethod] -> ShowS
show :: RootMethod -> String
$cshow :: RootMethod -> String
showsPrec :: Int -> RootMethod -> ShowS
$cshowsPrec :: Int -> RootMethod -> ShowS
Show,RootMethod
RootMethod -> RootMethod -> Bounded RootMethod
forall a. a -> a -> Bounded a
maxBound :: RootMethod
$cmaxBound :: RootMethod
minBound :: RootMethod
$cminBound :: RootMethod
Bounded)

-- | Nonlinear multidimensional root finding using algorithms that do not require 
-- any derivative information to be supplied by the user.
-- Any derivatives needed are approximated by finite differences.
root :: RootMethod
     -> Double                     -- ^ maximum residual
     -> Int                        -- ^ maximum number of iterations allowed
     -> ([Double] -> [Double])     -- ^ function to minimize
     -> [Double]                   -- ^ starting point
     -> ([Double], Matrix Double)  -- ^ solution vector and optimization path

root :: RootMethod
-> Double
-> Int
-> ([Double] -> [Double])
-> [Double]
-> ([Double], Matrix Double)
root method :: RootMethod
method epsabs :: Double
epsabs maxit :: Int
maxit fun :: [Double] -> [Double]
fun xinit :: [Double]
xinit = CInt
-> ([Double] -> [Double])
-> [Double]
-> Double
-> Int
-> ([Double], Matrix Double)
rootGen (Int -> CInt
fi (RootMethod -> Int
forall a. Enum a => a -> Int
fromEnum RootMethod
method)) [Double] -> [Double]
fun [Double]
xinit Double
epsabs Int
maxit

rootGen :: CInt
-> ([Double] -> [Double])
-> [Double]
-> Double
-> Int
-> ([Double], Matrix Double)
rootGen m :: CInt
m f :: [Double] -> [Double]
f xi :: [Double]
xi epsabs :: Double
epsabs maxit :: Int
maxit = IO ([Double], Matrix Double) -> ([Double], Matrix Double)
forall a. IO a -> a
unsafePerformIO (IO ([Double], Matrix Double) -> ([Double], Matrix Double))
-> IO ([Double], Matrix Double) -> ([Double], Matrix Double)
forall a b. (a -> b) -> a -> b
$ do
    let xiv :: Vector Double
xiv = [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList [Double]
xi
        n :: IndexOf Vector
n   = Vector Double -> IndexOf Vector
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
size Vector Double
xiv
    FunPtr TVV
fp <- TVV -> IO (FunPtr TVV)
mkVecVecfun ((Vector Double -> Vector Double) -> TVV
aux_vTov (IndexOf Vector -> Vector Double -> Vector Double
forall (c :: * -> *) t.
(Eq (IndexOf c), Container c t, Show (IndexOf c)) =>
IndexOf c -> c t -> c t
checkdim1 Int
IndexOf Vector
n (Vector Double -> Vector Double)
-> (Vector Double -> Vector Double)
-> Vector Double
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList ([Double] -> Vector Double)
-> (Vector Double -> [Double]) -> Vector Double -> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> [Double]
f ([Double] -> [Double])
-> (Vector Double -> [Double]) -> Vector Double -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
toList))
    Matrix Double
rawpath <- Vector Double
-> (((CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
     -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> IO (Matrix Double))
-> IO (Matrix Double)
forall a t b.
Storable a =>
Vector a -> (((CInt -> Ptr a -> t) -> t) -> IO b) -> IO b
vec Vector Double
xiv ((((CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
   -> CInt -> CInt -> Ptr Double -> IO CInt)
  -> IO (Matrix Double))
 -> IO (Matrix Double))
-> (((CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
     -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> IO (Matrix Double))
-> IO (Matrix Double)
forall a b. (a -> b) -> a -> b
$ \xiv' :: (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
-> CInt -> CInt -> Ptr Double -> IO CInt
xiv' ->
                   Int
-> Int
-> (CInt -> CInt -> Ptr Double -> IO CInt)
-> String
-> IO (Matrix Double)
forall a.
Storable a =>
Int
-> Int
-> (CInt -> CInt -> Ptr a -> IO CInt)
-> String
-> IO (Matrix a)
createMIO Int
maxit (2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+1)
                         (CInt
-> FunPtr TVV
-> Double
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
c_multiroot CInt
m FunPtr TVV
fp Double
epsabs (Int -> CInt
fi Int
maxit) (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
-> ((CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> CInt -> CInt -> Ptr Double -> IO CInt)
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
forall x y. x -> (x -> y) -> y
// (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
-> CInt -> CInt -> Ptr Double -> IO CInt
xiv')
                         "multiroot"
    let it :: Int
it = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Matrix Double
rawpath Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (Int
maxitInt -> Int -> Int
forall a. Num a => a -> a -> a
-1,0))
        path :: Matrix Double
path = Int -> Matrix Double -> Matrix Double
forall t. Element t => Int -> Matrix t -> Matrix t
takeRows Int
it Matrix Double
rawpath
        [sol :: [Double]
sol] = Matrix Double -> [[Double]]
forall t. Element t => Matrix t -> [[t]]
toLists (Matrix Double -> [[Double]]) -> Matrix Double -> [[Double]]
forall a b. (a -> b) -> a -> b
$ Int -> Matrix Double -> Matrix Double
forall t. Element t => Int -> Matrix t -> Matrix t
dropRows (Int
itInt -> Int -> Int
forall a. Num a => a -> a -> a
-1) Matrix Double
path
    FunPtr TVV -> IO ()
forall a. FunPtr a -> IO ()
freeHaskellFunPtr FunPtr TVV
fp
    ([Double], Matrix Double) -> IO ([Double], Matrix Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
n ([Double] -> [Double]) -> [Double] -> [Double]
forall a b. (a -> b) -> a -> b
$ Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop 1 [Double]
sol, Matrix Double
path)


foreign import ccall safe "multiroot"
    c_multiroot:: CInt -> FunPtr TVV -> Double -> CInt -> TVM

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

data RootMethodJ = HybridsJ
                 | HybridJ
                 | Newton
                 | GNewton
                deriving (Int -> RootMethodJ
RootMethodJ -> Int
RootMethodJ -> [RootMethodJ]
RootMethodJ -> RootMethodJ
RootMethodJ -> RootMethodJ -> [RootMethodJ]
RootMethodJ -> RootMethodJ -> RootMethodJ -> [RootMethodJ]
(RootMethodJ -> RootMethodJ)
-> (RootMethodJ -> RootMethodJ)
-> (Int -> RootMethodJ)
-> (RootMethodJ -> Int)
-> (RootMethodJ -> [RootMethodJ])
-> (RootMethodJ -> RootMethodJ -> [RootMethodJ])
-> (RootMethodJ -> RootMethodJ -> [RootMethodJ])
-> (RootMethodJ -> RootMethodJ -> RootMethodJ -> [RootMethodJ])
-> Enum RootMethodJ
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
enumFromThenTo :: RootMethodJ -> RootMethodJ -> RootMethodJ -> [RootMethodJ]
$cenumFromThenTo :: RootMethodJ -> RootMethodJ -> RootMethodJ -> [RootMethodJ]
enumFromTo :: RootMethodJ -> RootMethodJ -> [RootMethodJ]
$cenumFromTo :: RootMethodJ -> RootMethodJ -> [RootMethodJ]
enumFromThen :: RootMethodJ -> RootMethodJ -> [RootMethodJ]
$cenumFromThen :: RootMethodJ -> RootMethodJ -> [RootMethodJ]
enumFrom :: RootMethodJ -> [RootMethodJ]
$cenumFrom :: RootMethodJ -> [RootMethodJ]
fromEnum :: RootMethodJ -> Int
$cfromEnum :: RootMethodJ -> Int
toEnum :: Int -> RootMethodJ
$ctoEnum :: Int -> RootMethodJ
pred :: RootMethodJ -> RootMethodJ
$cpred :: RootMethodJ -> RootMethodJ
succ :: RootMethodJ -> RootMethodJ
$csucc :: RootMethodJ -> RootMethodJ
Enum,RootMethodJ -> RootMethodJ -> Bool
(RootMethodJ -> RootMethodJ -> Bool)
-> (RootMethodJ -> RootMethodJ -> Bool) -> Eq RootMethodJ
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: RootMethodJ -> RootMethodJ -> Bool
$c/= :: RootMethodJ -> RootMethodJ -> Bool
== :: RootMethodJ -> RootMethodJ -> Bool
$c== :: RootMethodJ -> RootMethodJ -> Bool
Eq,Int -> RootMethodJ -> ShowS
[RootMethodJ] -> ShowS
RootMethodJ -> String
(Int -> RootMethodJ -> ShowS)
-> (RootMethodJ -> String)
-> ([RootMethodJ] -> ShowS)
-> Show RootMethodJ
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [RootMethodJ] -> ShowS
$cshowList :: [RootMethodJ] -> ShowS
show :: RootMethodJ -> String
$cshow :: RootMethodJ -> String
showsPrec :: Int -> RootMethodJ -> ShowS
$cshowsPrec :: Int -> RootMethodJ -> ShowS
Show,RootMethodJ
RootMethodJ -> RootMethodJ -> Bounded RootMethodJ
forall a. a -> a -> Bounded a
maxBound :: RootMethodJ
$cmaxBound :: RootMethodJ
minBound :: RootMethodJ
$cminBound :: RootMethodJ
Bounded)

-- | Nonlinear multidimensional root finding using both the function and its derivatives.
rootJ :: RootMethodJ
      -> Double                     -- ^ maximum residual
      -> Int                        -- ^ maximum number of iterations allowed
      -> ([Double] -> [Double])     -- ^ function to minimize
      -> ([Double] -> [[Double]])   -- ^ Jacobian
      -> [Double]                   -- ^ starting point
      -> ([Double], Matrix Double)  -- ^ solution vector and optimization path

rootJ :: RootMethodJ
-> Double
-> Int
-> ([Double] -> [Double])
-> ([Double] -> [[Double]])
-> [Double]
-> ([Double], Matrix Double)
rootJ method :: RootMethodJ
method epsabs :: Double
epsabs maxit :: Int
maxit fun :: [Double] -> [Double]
fun jac :: [Double] -> [[Double]]
jac xinit :: [Double]
xinit = CInt
-> ([Double] -> [Double])
-> ([Double] -> [[Double]])
-> [Double]
-> Double
-> Int
-> ([Double], Matrix Double)
rootJGen (Int -> CInt
fi (RootMethodJ -> Int
forall a. Enum a => a -> Int
fromEnum RootMethodJ
method)) [Double] -> [Double]
fun [Double] -> [[Double]]
jac [Double]
xinit Double
epsabs Int
maxit

rootJGen :: CInt
-> ([Double] -> [Double])
-> ([Double] -> [[Double]])
-> [Double]
-> Double
-> Int
-> ([Double], Matrix Double)
rootJGen m :: CInt
m f :: [Double] -> [Double]
f jac :: [Double] -> [[Double]]
jac xi :: [Double]
xi epsabs :: Double
epsabs maxit :: Int
maxit = IO ([Double], Matrix Double) -> ([Double], Matrix Double)
forall a. IO a -> a
unsafePerformIO (IO ([Double], Matrix Double) -> ([Double], Matrix Double))
-> IO ([Double], Matrix Double) -> ([Double], Matrix Double)
forall a b. (a -> b) -> a -> b
$ do
    let xiv :: Vector Double
xiv = [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList [Double]
xi
        n :: IndexOf Vector
n   = Vector Double -> IndexOf Vector
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
size Vector Double
xiv
    FunPtr TVV
fp <- TVV -> IO (FunPtr TVV)
mkVecVecfun ((Vector Double -> Vector Double) -> TVV
aux_vTov (IndexOf Vector -> Vector Double -> Vector Double
forall (c :: * -> *) t.
(Eq (IndexOf c), Container c t, Show (IndexOf c)) =>
IndexOf c -> c t -> c t
checkdim1 Int
IndexOf Vector
n (Vector Double -> Vector Double)
-> (Vector Double -> Vector Double)
-> Vector Double
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList ([Double] -> Vector Double)
-> (Vector Double -> [Double]) -> Vector Double -> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> [Double]
f ([Double] -> [Double])
-> (Vector Double -> [Double]) -> Vector Double -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
toList))
    FunPtr
  (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
jp <- (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
-> IO
     (FunPtr
        (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt))
mkVecMatfun ((Vector Double -> Matrix Double)
-> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt
aux_vTom (Int -> Matrix Double -> Matrix Double
forall t. Int -> Matrix t -> Matrix t
checkdim2 Int
n (Matrix Double -> Matrix Double)
-> (Vector Double -> Matrix Double)
-> Vector Double
-> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [[Double]] -> Matrix Double
forall t. Element t => [[t]] -> Matrix t
fromLists ([[Double]] -> Matrix Double)
-> (Vector Double -> [[Double]]) -> Vector Double -> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> [[Double]]
jac ([Double] -> [[Double]])
-> (Vector Double -> [Double]) -> Vector Double -> [[Double]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
toList))
    Matrix Double
rawpath <- Vector Double
-> (((CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
     -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> IO (Matrix Double))
-> IO (Matrix Double)
forall a t b.
Storable a =>
Vector a -> (((CInt -> Ptr a -> t) -> t) -> IO b) -> IO b
vec Vector Double
xiv ((((CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
   -> CInt -> CInt -> Ptr Double -> IO CInt)
  -> IO (Matrix Double))
 -> IO (Matrix Double))
-> (((CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
     -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> IO (Matrix Double))
-> IO (Matrix Double)
forall a b. (a -> b) -> a -> b
$ \xiv' :: (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
-> CInt -> CInt -> Ptr Double -> IO CInt
xiv' ->
                   Int
-> Int
-> (CInt -> CInt -> Ptr Double -> IO CInt)
-> String
-> IO (Matrix Double)
forall a.
Storable a =>
Int
-> Int
-> (CInt -> CInt -> Ptr a -> IO CInt)
-> String
-> IO (Matrix a)
createMIO Int
maxit (2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+1)
                         (CInt
-> FunPtr TVV
-> FunPtr
     (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
-> Double
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
c_multirootj CInt
m FunPtr TVV
fp FunPtr
  (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
jp Double
epsabs (Int -> CInt
fi Int
maxit) (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
-> ((CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> CInt -> CInt -> Ptr Double -> IO CInt)
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
forall x y. x -> (x -> y) -> y
// (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
-> CInt -> CInt -> Ptr Double -> IO CInt
xiv')
                         "multiroot"
    let it :: Int
it = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Matrix Double
rawpath Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (Int
maxitInt -> Int -> Int
forall a. Num a => a -> a -> a
-1,0))
        path :: Matrix Double
path = Int -> Matrix Double -> Matrix Double
forall t. Element t => Int -> Matrix t -> Matrix t
takeRows Int
it Matrix Double
rawpath
        [sol :: [Double]
sol] = Matrix Double -> [[Double]]
forall t. Element t => Matrix t -> [[t]]
toLists (Matrix Double -> [[Double]]) -> Matrix Double -> [[Double]]
forall a b. (a -> b) -> a -> b
$ Int -> Matrix Double -> Matrix Double
forall t. Element t => Int -> Matrix t -> Matrix t
dropRows (Int
itInt -> Int -> Int
forall a. Num a => a -> a -> a
-1) Matrix Double
path
    FunPtr TVV -> IO ()
forall a. FunPtr a -> IO ()
freeHaskellFunPtr FunPtr TVV
fp
    FunPtr
  (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
-> IO ()
forall a. FunPtr a -> IO ()
freeHaskellFunPtr FunPtr
  (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
jp
    ([Double], Matrix Double) -> IO ([Double], Matrix Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
n ([Double] -> [Double]) -> [Double] -> [Double]
forall a b. (a -> b) -> a -> b
$ Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop 1 [Double]
sol, Matrix Double
path)

foreign import ccall safe "multirootj"
    c_multirootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM

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

checkdim1 :: IndexOf c -> c t -> c t
checkdim1 n :: IndexOf c
n v :: c t
v
    | c t -> IndexOf c
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
size c t
v IndexOf c -> IndexOf c -> Bool
forall a. Eq a => a -> a -> Bool
== IndexOf c
n = c t
v
    | Bool
otherwise = String -> c t
forall a. HasCallStack => String -> a
error (String -> c t) -> String -> c t
forall a b. (a -> b) -> a -> b
$ "Error: "String -> ShowS
forall a. [a] -> [a] -> [a]
++ IndexOf c -> String
forall a. Show a => a -> String
show IndexOf c
n
                        String -> ShowS
forall a. [a] -> [a] -> [a]
++ " components expected in the result of the function supplied to root"

checkdim2 :: Int -> Matrix t -> Matrix t
checkdim2 n :: Int
n m :: Matrix t
m
    | Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n Bool -> Bool -> Bool
&& Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n = Matrix t
m
    | Bool
otherwise = String -> Matrix t
forall a. HasCallStack => String -> a
error (String -> Matrix t) -> String -> Matrix t
forall a b. (a -> b) -> a -> b
$ "Error: "String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show Int
n String -> ShowS
forall a. [a] -> [a] -> [a]
++ "x" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show Int
n
                        String -> ShowS
forall a. [a] -> [a] -> [a]
++ " Jacobian expected in rootJ"