diff options
author | Alex Lang <me@alang.ca> | 2013-05-08 10:36:47 +0900 |
---|---|---|
committer | Alex Lang <me@alang.ca> | 2013-05-08 10:36:47 +0900 |
commit | edb20f3c9993d9c9c349f6a001317592949154f1 (patch) | |
tree | 3ac07f1d4391529808c98d85aba9f77859499217 | |
parent | 17c56e5f9563245eb3aaa042c843745fce096554 (diff) |
Implemented uniRoot for one-dimensional root-finding without derivatives
-rw-r--r-- | lib/Numeric/GSL/Internal.hs | 3 | ||||
-rw-r--r-- | lib/Numeric/GSL/Root.hs | 44 | ||||
-rw-r--r-- | lib/Numeric/GSL/gsl-aux.c | 65 |
3 files changed, 101 insertions, 11 deletions
diff --git a/lib/Numeric/GSL/Internal.hs b/lib/Numeric/GSL/Internal.hs index 84417ce..69a9750 100644 --- a/lib/Numeric/GSL/Internal.hs +++ b/lib/Numeric/GSL/Internal.hs | |||
@@ -36,6 +36,9 @@ foreign import ccall safe "wrapper" | |||
36 | foreign import ccall safe "wrapper" | 36 | foreign import ccall safe "wrapper" |
37 | mkDoubleVecVecfun :: (Double -> TVV) -> IO (FunPtr (Double -> TVV)) | 37 | mkDoubleVecVecfun :: (Double -> TVV) -> IO (FunPtr (Double -> TVV)) |
38 | 38 | ||
39 | foreign import ccall safe "wrapper" | ||
40 | mkDoublefun :: (Double -> Double) -> IO (FunPtr (Double -> Double)) | ||
41 | |||
39 | aux_vTov :: (Vector Double -> Vector Double) -> TVV | 42 | aux_vTov :: (Vector Double -> Vector Double) -> TVV |
40 | aux_vTov f n p nr r = g where | 43 | aux_vTov f n p nr r = g where |
41 | v = f x | 44 | v = f x |
diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs index cd2982a..25ec5f5 100644 --- a/lib/Numeric/GSL/Root.hs +++ b/lib/Numeric/GSL/Root.hs | |||
@@ -45,6 +45,7 @@ main = do | |||
45 | ----------------------------------------------------------------------------- | 45 | ----------------------------------------------------------------------------- |
46 | 46 | ||
47 | module Numeric.GSL.Root ( | 47 | module Numeric.GSL.Root ( |
48 | uniRoot, UniRootMethod(..), | ||
48 | root, RootMethod(..), | 49 | root, RootMethod(..), |
49 | rootJ, RootMethodJ(..), | 50 | rootJ, RootMethodJ(..), |
50 | ) where | 51 | ) where |
@@ -58,6 +59,36 @@ import System.IO.Unsafe(unsafePerformIO) | |||
58 | 59 | ||
59 | ------------------------------------------------------------------------- | 60 | ------------------------------------------------------------------------- |
60 | 61 | ||
62 | data UniRootMethod = Bisection | ||
63 | | FalsePos | ||
64 | | Brent | ||
65 | deriving (Enum, Eq, Show, Bounded) | ||
66 | |||
67 | uniRoot :: UniRootMethod | ||
68 | -> Double | ||
69 | -> Int | ||
70 | -> (Double -> Double) | ||
71 | -> Double | ||
72 | -> Double | ||
73 | -> (Double, Matrix Double) | ||
74 | uniRoot method epsrel maxit fun xl xu = uniRootGen (fi (fromEnum method)) fun xl xu epsrel maxit | ||
75 | |||
76 | uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do | ||
77 | fp <- mkDoublefun f | ||
78 | rawpath <- createMIO maxit 4 | ||
79 | (c_root m fp epsrel (fi maxit) xl xu) | ||
80 | "root" | ||
81 | let it = round (rawpath @@> (maxit-1,0)) | ||
82 | path = takeRows it rawpath | ||
83 | [sol] = toLists $ dropRows (it-1) path | ||
84 | freeHaskellFunPtr fp | ||
85 | return (sol !! 1, path) | ||
86 | |||
87 | foreign import ccall safe "root" | ||
88 | c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM | ||
89 | |||
90 | ------------------------------------------------------------------------- | ||
91 | |||
61 | data RootMethod = Hybrids | 92 | data RootMethod = Hybrids |
62 | | Hybrid | 93 | | Hybrid |
63 | | DNewton | 94 | | DNewton |
@@ -82,8 +113,8 @@ rootGen m f xi epsabs maxit = unsafePerformIO $ do | |||
82 | fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) | 113 | fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) |
83 | rawpath <- vec xiv $ \xiv' -> | 114 | rawpath <- vec xiv $ \xiv' -> |
84 | createMIO maxit (2*n+1) | 115 | createMIO maxit (2*n+1) |
85 | (c_root m fp epsabs (fi maxit) // xiv') | 116 | (c_multiroot m fp epsabs (fi maxit) // xiv') |
86 | "root" | 117 | "multiroot" |
87 | let it = round (rawpath @@> (maxit-1,0)) | 118 | let it = round (rawpath @@> (maxit-1,0)) |
88 | path = takeRows it rawpath | 119 | path = takeRows it rawpath |
89 | [sol] = toLists $ dropRows (it-1) path | 120 | [sol] = toLists $ dropRows (it-1) path |
@@ -91,8 +122,8 @@ rootGen m f xi epsabs maxit = unsafePerformIO $ do | |||
91 | return (take n $ drop 1 sol, path) | 122 | return (take n $ drop 1 sol, path) |
92 | 123 | ||
93 | 124 | ||
94 | foreign import ccall safe "root" | 125 | foreign import ccall safe "multiroot" |
95 | c_root:: CInt -> FunPtr TVV -> Double -> CInt -> TVM | 126 | c_multiroot:: CInt -> FunPtr TVV -> Double -> CInt -> TVM |
96 | 127 | ||
97 | ------------------------------------------------------------------------- | 128 | ------------------------------------------------------------------------- |
98 | 129 | ||
@@ -121,7 +152,7 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do | |||
121 | rawpath <- vec xiv $ \xiv' -> | 152 | rawpath <- vec xiv $ \xiv' -> |
122 | createMIO maxit (2*n+1) | 153 | createMIO maxit (2*n+1) |
123 | (c_rootj m fp jp epsabs (fi maxit) // xiv') | 154 | (c_rootj m fp jp epsabs (fi maxit) // xiv') |
124 | "root" | 155 | "multiroot" |
125 | let it = round (rawpath @@> (maxit-1,0)) | 156 | let it = round (rawpath @@> (maxit-1,0)) |
126 | path = takeRows it rawpath | 157 | path = takeRows it rawpath |
127 | [sol] = toLists $ dropRows (it-1) path | 158 | [sol] = toLists $ dropRows (it-1) path |
@@ -129,8 +160,7 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do | |||
129 | freeHaskellFunPtr jp | 160 | freeHaskellFunPtr jp |
130 | return (take n $ drop 1 sol, path) | 161 | return (take n $ drop 1 sol, path) |
131 | 162 | ||
132 | 163 | foreign import ccall safe "multirootj" | |
133 | foreign import ccall safe "rootj" | ||
134 | c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM | 164 | c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM |
135 | 165 | ||
136 | ------------------------------------------------------- | 166 | ------------------------------------------------------- |
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c index 61b6a54..27ffbdc 100644 --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c | |||
@@ -32,6 +32,7 @@ | |||
32 | #include <gsl/gsl_complex_math.h> | 32 | #include <gsl/gsl_complex_math.h> |
33 | #include <gsl/gsl_rng.h> | 33 | #include <gsl/gsl_rng.h> |
34 | #include <gsl/gsl_randist.h> | 34 | #include <gsl/gsl_randist.h> |
35 | #include <gsl/gsl_roots.h> | ||
35 | #include <gsl/gsl_multifit_nlin.h> | 36 | #include <gsl/gsl_multifit_nlin.h> |
36 | #include <string.h> | 37 | #include <string.h> |
37 | #include <stdio.h> | 38 | #include <stdio.h> |
@@ -1047,9 +1048,65 @@ int minimizeD(int method, double f(int, double*), int df(int, double*, int, doub | |||
1047 | 1048 | ||
1048 | //--------------------------------------------------------------- | 1049 | //--------------------------------------------------------------- |
1049 | 1050 | ||
1051 | double only_f_aux_root(double x, void *pars) { | ||
1052 | double (*f)(double) = (double (*)(double)) pars; | ||
1053 | return f(x); | ||
1054 | } | ||
1055 | |||
1056 | int root(int method, double f(double), | ||
1057 | double epsrel, int maxit, | ||
1058 | double xl, double xu, RMAT(sol)) { | ||
1059 | REQUIRES(solr == maxit && solc == 4,BAD_SIZE); | ||
1060 | DEBUGMSG("root_only_f"); | ||
1061 | gsl_function my_func; | ||
1062 | // extract function from pars | ||
1063 | my_func.function = only_f_aux_root; | ||
1064 | my_func.params = f; | ||
1065 | size_t iter = 0; | ||
1066 | int status; | ||
1067 | const gsl_root_fsolver_type *T; | ||
1068 | gsl_root_fsolver *s; | ||
1069 | // Starting point | ||
1070 | switch(method) { | ||
1071 | case 0 : {T = gsl_root_fsolver_bisection; printf("7\n"); break; } | ||
1072 | case 1 : {T = gsl_root_fsolver_falsepos; break; } | ||
1073 | case 2 : {T = gsl_root_fsolver_brent; break; } | ||
1074 | default: ERROR(BAD_CODE); | ||
1075 | } | ||
1076 | s = gsl_root_fsolver_alloc (T); | ||
1077 | gsl_root_fsolver_set (s, &my_func, xl, xu); | ||
1078 | do { | ||
1079 | double best, current_lo, current_hi; | ||
1080 | status = gsl_root_fsolver_iterate (s); | ||
1081 | best = gsl_root_fsolver_root (s); | ||
1082 | current_lo = gsl_root_fsolver_x_lower (s); | ||
1083 | current_hi = gsl_root_fsolver_x_upper (s); | ||
1084 | solp[iter*solc] = iter + 1; | ||
1085 | solp[iter*solc+1] = best; | ||
1086 | solp[iter*solc+2] = current_lo; | ||
1087 | solp[iter*solc+3] = current_hi; | ||
1088 | iter++; | ||
1089 | if (status) /* check if solver is stuck */ | ||
1090 | break; | ||
1091 | |||
1092 | status = | ||
1093 | gsl_root_test_interval (current_lo, current_hi, 0, epsrel); | ||
1094 | } | ||
1095 | while (status == GSL_CONTINUE && iter < maxit); | ||
1096 | int i; | ||
1097 | for (i=iter; i<solr; i++) { | ||
1098 | solp[i*solc+0] = iter; | ||
1099 | solp[i*solc+1]=0.; | ||
1100 | solp[i*solc+2]=0.; | ||
1101 | solp[i*solc+3]=0.; | ||
1102 | } | ||
1103 | gsl_root_fsolver_free(s); | ||
1104 | OK | ||
1105 | } | ||
1106 | |||
1050 | typedef void TrawfunV(int, double*, int, double*); | 1107 | typedef void TrawfunV(int, double*, int, double*); |
1051 | 1108 | ||
1052 | int only_f_aux_root(const gsl_vector*x, void *pars, gsl_vector*y) { | 1109 | int only_f_aux_multiroot(const gsl_vector*x, void *pars, gsl_vector*y) { |
1053 | TrawfunV * f = (TrawfunV*) pars; | 1110 | TrawfunV * f = (TrawfunV*) pars; |
1054 | double* p = (double*)calloc(x->size,sizeof(double)); | 1111 | double* p = (double*)calloc(x->size,sizeof(double)); |
1055 | double* q = (double*)calloc(y->size,sizeof(double)); | 1112 | double* q = (double*)calloc(y->size,sizeof(double)); |
@@ -1066,14 +1123,14 @@ int only_f_aux_root(const gsl_vector*x, void *pars, gsl_vector*y) { | |||
1066 | return 0; //hmmm | 1123 | return 0; //hmmm |
1067 | } | 1124 | } |
1068 | 1125 | ||
1069 | int root(int method, void f(int, double*, int, double*), | 1126 | int multiroot(int method, void f(int, double*, int, double*), |
1070 | double epsabs, int maxit, | 1127 | double epsabs, int maxit, |
1071 | KRVEC(xi), RMAT(sol)) { | 1128 | KRVEC(xi), RMAT(sol)) { |
1072 | REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE); | 1129 | REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE); |
1073 | DEBUGMSG("root_only_f"); | 1130 | DEBUGMSG("root_only_f"); |
1074 | gsl_multiroot_function my_func; | 1131 | gsl_multiroot_function my_func; |
1075 | // extract function from pars | 1132 | // extract function from pars |
1076 | my_func.f = only_f_aux_root; | 1133 | my_func.f = only_f_aux_multiroot; |
1077 | my_func.n = xin; | 1134 | my_func.n = xin; |
1078 | my_func.params = f; | 1135 | my_func.params = f; |
1079 | size_t iter = 0; | 1136 | size_t iter = 0; |
@@ -1175,7 +1232,7 @@ int fjf_aux(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) { | |||
1175 | return 0; | 1232 | return 0; |
1176 | } | 1233 | } |
1177 | 1234 | ||
1178 | int rootj(int method, int f(int, double*, int, double*), | 1235 | int multirootj(int method, int f(int, double*, int, double*), |
1179 | int jac(int, double*, int, int, double*), | 1236 | int jac(int, double*, int, int, double*), |
1180 | double epsabs, int maxit, | 1237 | double epsabs, int maxit, |
1181 | KRVEC(xi), RMAT(sol)) { | 1238 | KRVEC(xi), RMAT(sol)) { |