{-# 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
$csucc :: UniRootMethod -> UniRootMethod
succ :: UniRootMethod -> UniRootMethod
$cpred :: UniRootMethod -> UniRootMethod
pred :: UniRootMethod -> UniRootMethod
$ctoEnum :: Int -> UniRootMethod
toEnum :: Int -> UniRootMethod
$cfromEnum :: UniRootMethod -> Int
fromEnum :: UniRootMethod -> Int
$cenumFrom :: UniRootMethod -> [UniRootMethod]
enumFrom :: UniRootMethod -> [UniRootMethod]
$cenumFromThen :: UniRootMethod -> UniRootMethod -> [UniRootMethod]
enumFromThen :: UniRootMethod -> UniRootMethod -> [UniRootMethod]
$cenumFromTo :: UniRootMethod -> UniRootMethod -> [UniRootMethod]
enumFromTo :: UniRootMethod -> UniRootMethod -> [UniRootMethod]
$cenumFromThenTo :: UniRootMethod -> UniRootMethod -> UniRootMethod -> [UniRootMethod]
enumFromThenTo :: UniRootMethod -> UniRootMethod -> UniRootMethod -> [UniRootMethod]
Enum, UniRootMethod -> UniRootMethod -> Bool
(UniRootMethod -> UniRootMethod -> Bool)
-> (UniRootMethod -> UniRootMethod -> Bool) -> Eq UniRootMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: UniRootMethod -> UniRootMethod -> Bool
== :: UniRootMethod -> UniRootMethod -> Bool
$c/= :: UniRootMethod -> UniRootMethod -> Bool
/= :: 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
$cshowsPrec :: Int -> UniRootMethod -> ShowS
showsPrec :: Int -> UniRootMethod -> ShowS
$cshow :: UniRootMethod -> String
show :: UniRootMethod -> String
$cshowList :: [UniRootMethod] -> ShowS
showList :: [UniRootMethod] -> ShowS
Show, UniRootMethod
UniRootMethod -> UniRootMethod -> Bounded UniRootMethod
forall a. a -> a -> Bounded a
$cminBound :: UniRootMethod
minBound :: UniRootMethod
$cmaxBound :: UniRootMethod
maxBound :: 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 UniRootMethod
method Double
epsrel Int
maxit Double -> Double
fun Double
xl 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 CInt
m Double -> Double
f Double
xl Double
xu Double
epsrel 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 Int
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)
                         String
"root"
    let it :: Int
it = Double -> Int
forall b. Integral b => Double -> b
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
-Int
1,Int
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
        [[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
-Int
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 a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return ([Double]
sol [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
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
$csucc :: UniRootMethodJ -> UniRootMethodJ
succ :: UniRootMethodJ -> UniRootMethodJ
$cpred :: UniRootMethodJ -> UniRootMethodJ
pred :: UniRootMethodJ -> UniRootMethodJ
$ctoEnum :: Int -> UniRootMethodJ
toEnum :: Int -> UniRootMethodJ
$cfromEnum :: UniRootMethodJ -> Int
fromEnum :: UniRootMethodJ -> Int
$cenumFrom :: UniRootMethodJ -> [UniRootMethodJ]
enumFrom :: UniRootMethodJ -> [UniRootMethodJ]
$cenumFromThen :: UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
enumFromThen :: UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
$cenumFromTo :: UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
enumFromTo :: UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
$cenumFromThenTo :: UniRootMethodJ
-> UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
enumFromThenTo :: UniRootMethodJ
-> UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
Enum, UniRootMethodJ -> UniRootMethodJ -> Bool
(UniRootMethodJ -> UniRootMethodJ -> Bool)
-> (UniRootMethodJ -> UniRootMethodJ -> Bool) -> Eq UniRootMethodJ
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: UniRootMethodJ -> UniRootMethodJ -> Bool
== :: UniRootMethodJ -> UniRootMethodJ -> Bool
$c/= :: UniRootMethodJ -> UniRootMethodJ -> Bool
/= :: 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
$cshowsPrec :: Int -> UniRootMethodJ -> ShowS
showsPrec :: Int -> UniRootMethodJ -> ShowS
$cshow :: UniRootMethodJ -> String
show :: UniRootMethodJ -> String
$cshowList :: [UniRootMethodJ] -> ShowS
showList :: [UniRootMethodJ] -> ShowS
Show, UniRootMethodJ
UniRootMethodJ -> UniRootMethodJ -> Bounded UniRootMethodJ
forall a. a -> a -> Bounded a
$cminBound :: UniRootMethodJ
minBound :: UniRootMethodJ
$cmaxBound :: UniRootMethodJ
maxBound :: 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 UniRootMethodJ
method Double
epsrel Int
maxit Double -> Double
fun Double -> Double
dfun 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 CInt
m Double -> Double
f Double -> Double
df Double
x Double
epsrel 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 Int
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)
                         String
"rootj"
    let it :: Int
it = Double -> Int
forall b. Integral b => Double -> b
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
-Int
1,Int
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
        [[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
-Int
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 a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return ([Double]
sol [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
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
$csucc :: RootMethod -> RootMethod
succ :: RootMethod -> RootMethod
$cpred :: RootMethod -> RootMethod
pred :: RootMethod -> RootMethod
$ctoEnum :: Int -> RootMethod
toEnum :: Int -> RootMethod
$cfromEnum :: RootMethod -> Int
fromEnum :: RootMethod -> Int
$cenumFrom :: RootMethod -> [RootMethod]
enumFrom :: RootMethod -> [RootMethod]
$cenumFromThen :: RootMethod -> RootMethod -> [RootMethod]
enumFromThen :: RootMethod -> RootMethod -> [RootMethod]
$cenumFromTo :: RootMethod -> RootMethod -> [RootMethod]
enumFromTo :: RootMethod -> RootMethod -> [RootMethod]
$cenumFromThenTo :: RootMethod -> RootMethod -> RootMethod -> [RootMethod]
enumFromThenTo :: RootMethod -> RootMethod -> RootMethod -> [RootMethod]
Enum,RootMethod -> RootMethod -> Bool
(RootMethod -> RootMethod -> Bool)
-> (RootMethod -> RootMethod -> Bool) -> Eq RootMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RootMethod -> RootMethod -> Bool
== :: RootMethod -> RootMethod -> Bool
$c/= :: RootMethod -> RootMethod -> Bool
/= :: 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
$cshowsPrec :: Int -> RootMethod -> ShowS
showsPrec :: Int -> RootMethod -> ShowS
$cshow :: RootMethod -> String
show :: RootMethod -> String
$cshowList :: [RootMethod] -> ShowS
showList :: [RootMethod] -> ShowS
Show,RootMethod
RootMethod -> RootMethod -> Bounded RootMethod
forall a. a -> a -> Bounded a
$cminBound :: RootMethod
minBound :: RootMethod
$cmaxBound :: RootMethod
maxBound :: 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 RootMethod
method Double
epsabs Int
maxit [Double] -> [Double]
fun [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 CInt
m [Double] -> [Double]
f [Double]
xi Double
epsabs 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
$ \(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 (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
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')
                         String
"multiroot"
    let it :: Int
it = Double -> Int
forall b. Integral b => Double -> b
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
-Int
1,Int
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
        [[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
-Int
1) Matrix Double
path
    FunPtr TVV -> IO ()
forall a. FunPtr a -> IO ()
freeHaskellFunPtr FunPtr TVV
fp
    ([Double], Matrix Double) -> IO ([Double], Matrix Double)
forall a. a -> IO a
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 Int
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
$csucc :: RootMethodJ -> RootMethodJ
succ :: RootMethodJ -> RootMethodJ
$cpred :: RootMethodJ -> RootMethodJ
pred :: RootMethodJ -> RootMethodJ
$ctoEnum :: Int -> RootMethodJ
toEnum :: Int -> RootMethodJ
$cfromEnum :: RootMethodJ -> Int
fromEnum :: RootMethodJ -> Int
$cenumFrom :: RootMethodJ -> [RootMethodJ]
enumFrom :: RootMethodJ -> [RootMethodJ]
$cenumFromThen :: RootMethodJ -> RootMethodJ -> [RootMethodJ]
enumFromThen :: RootMethodJ -> RootMethodJ -> [RootMethodJ]
$cenumFromTo :: RootMethodJ -> RootMethodJ -> [RootMethodJ]
enumFromTo :: RootMethodJ -> RootMethodJ -> [RootMethodJ]
$cenumFromThenTo :: RootMethodJ -> RootMethodJ -> RootMethodJ -> [RootMethodJ]
enumFromThenTo :: RootMethodJ -> RootMethodJ -> RootMethodJ -> [RootMethodJ]
Enum,RootMethodJ -> RootMethodJ -> Bool
(RootMethodJ -> RootMethodJ -> Bool)
-> (RootMethodJ -> RootMethodJ -> Bool) -> Eq RootMethodJ
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RootMethodJ -> RootMethodJ -> Bool
== :: RootMethodJ -> RootMethodJ -> Bool
$c/= :: RootMethodJ -> RootMethodJ -> Bool
/= :: 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
$cshowsPrec :: Int -> RootMethodJ -> ShowS
showsPrec :: Int -> RootMethodJ -> ShowS
$cshow :: RootMethodJ -> String
show :: RootMethodJ -> String
$cshowList :: [RootMethodJ] -> ShowS
showList :: [RootMethodJ] -> ShowS
Show,RootMethodJ
RootMethodJ -> RootMethodJ -> Bounded RootMethodJ
forall a. a -> a -> Bounded a
$cminBound :: RootMethodJ
minBound :: RootMethodJ
$cmaxBound :: RootMethodJ
maxBound :: 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 RootMethodJ
method Double
epsabs Int
maxit [Double] -> [Double]
fun [Double] -> [[Double]]
jac [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 CInt
m [Double] -> [Double]
f [Double] -> [[Double]]
jac [Double]
xi Double
epsabs 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
$ \(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 (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
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')
                         String
"multiroot"
    let it :: Int
it = Double -> Int
forall b. Integral b => Double -> b
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
-Int
1,Int
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
        [[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
-Int
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 a. a -> IO a
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 Int
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 IndexOf c
n 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
$ String
"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]
++ String
" components expected in the result of the function supplied to root"

checkdim2 :: Int -> Matrix t -> Matrix t
checkdim2 Int
n 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
$ String
"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]
++ String
"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]
++ String
" Jacobian expected in rootJ"