{-# LANGUAGE FlexibleContexts #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
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 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 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 (m :: * -> *) a. Monad m => a -> m a
return ([Double]
sol [Double] -> Int -> Double
forall a. [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
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 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 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 (m :: * -> *) a. Monad m => a -> m a
return ([Double]
sol [Double] -> Int -> Double
forall a. [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
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)
root :: RootMethod
-> Double
-> Int
-> ([Double] -> [Double])
-> [Double]
-> ([Double], Matrix Double)
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 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 (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
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)
rootJ :: RootMethodJ
-> Double
-> Int
-> ([Double] -> [Double])
-> ([Double] -> [[Double]])
-> [Double]
-> ([Double], Matrix Double)
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 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 (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"