diff options
-rw-r--r-- | lib/Numeric/GSL/Minimization.hs | 34 | ||||
-rw-r--r-- | lib/Numeric/GSL/gsl-aux.c | 54 |
2 files changed, 88 insertions, 0 deletions
diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs index af85135..7ccca81 100644 --- a/lib/Numeric/GSL/Minimization.hs +++ b/lib/Numeric/GSL/Minimization.hs | |||
@@ -53,6 +53,7 @@ The nmsimplex2 version is a new O(N) implementation of the earlier O(N^2) nmsimp | |||
53 | module Numeric.GSL.Minimization ( | 53 | module Numeric.GSL.Minimization ( |
54 | minimize, minimizeV, MinimizeMethod(..), | 54 | minimize, minimizeV, MinimizeMethod(..), |
55 | minimizeD, minimizeVD, MinimizeMethodD(..), | 55 | minimizeD, minimizeVD, MinimizeMethodD(..), |
56 | uniMinimize, UniMinimizeMethod(..), | ||
56 | 57 | ||
57 | minimizeNMSimplex, | 58 | minimizeNMSimplex, |
58 | minimizeConjugateGradient, | 59 | minimizeConjugateGradient, |
@@ -81,6 +82,39 @@ minimizeVectorBFGS2 step tol eps maxit f g xi = minimizeD VectorBFGS2 eps maxit | |||
81 | 82 | ||
82 | ------------------------------------------------------------------------- | 83 | ------------------------------------------------------------------------- |
83 | 84 | ||
85 | data UniMinimizeMethod = GoldenSection | ||
86 | | BrentMini | ||
87 | | QuadGolden | ||
88 | deriving (Enum, Eq, Show, Bounded) | ||
89 | |||
90 | -- | Onedimensional minimization. | ||
91 | |||
92 | uniMinimize :: UniMinimizeMethod -- ^ The method used. | ||
93 | -> Double -- ^ desired precision of the solution | ||
94 | -> Int -- ^ maximum number of iterations allowed | ||
95 | -> (Double -> Double) -- ^ function to minimize | ||
96 | -> Double -- ^ guess for the location of the minimum | ||
97 | -> Double -- ^ lower bound of search interval | ||
98 | -> Double -- ^ upper bound of search interval | ||
99 | -> (Double, Matrix Double) -- ^ solution and optimization path | ||
100 | |||
101 | uniMinimize method epsrel maxit fun xmin xl xu = uniMinimizeGen (fi (fromEnum method)) fun xmin xl xu epsrel maxit | ||
102 | |||
103 | uniMinimizeGen m f xmin xl xu epsrel maxit = unsafePerformIO $ do | ||
104 | fp <- mkDoublefun f | ||
105 | rawpath <- createMIO maxit 4 | ||
106 | (c_uniMinize m fp epsrel (fi maxit) xmin xl xu) | ||
107 | "uniMinimize" | ||
108 | let it = round (rawpath @@> (maxit-1,0)) | ||
109 | path = takeRows it rawpath | ||
110 | [sol] = toLists $ dropRows (it-1) path | ||
111 | freeHaskellFunPtr fp | ||
112 | return (sol !! 1, path) | ||
113 | |||
114 | |||
115 | foreign import ccall safe "uniMinimize" | ||
116 | c_uniMinize:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> Double -> TM | ||
117 | |||
84 | data MinimizeMethod = NMSimplex | 118 | data MinimizeMethod = NMSimplex |
85 | | NMSimplex2 | 119 | | NMSimplex2 |
86 | deriving (Enum,Eq,Show,Bounded) | 120 | deriving (Enum,Eq,Show,Bounded) |
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c index e727c91..58bccb3 100644 --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c | |||
@@ -29,6 +29,7 @@ | |||
29 | #include <gsl/gsl_poly.h> | 29 | #include <gsl/gsl_poly.h> |
30 | #include <gsl/gsl_multimin.h> | 30 | #include <gsl/gsl_multimin.h> |
31 | #include <gsl/gsl_multiroots.h> | 31 | #include <gsl/gsl_multiroots.h> |
32 | #include <gsl/gsl_min.h> | ||
32 | #include <gsl/gsl_complex_math.h> | 33 | #include <gsl/gsl_complex_math.h> |
33 | #include <gsl/gsl_rng.h> | 34 | #include <gsl/gsl_rng.h> |
34 | #include <gsl/gsl_randist.h> | 35 | #include <gsl/gsl_randist.h> |
@@ -892,6 +893,59 @@ double only_f_aux_min(const gsl_vector*x, void *pars) { | |||
892 | return res; | 893 | return res; |
893 | } | 894 | } |
894 | 895 | ||
896 | double only_f_aux_root(double x, void *pars); | ||
897 | int uniMinimize(int method, double f(double), | ||
898 | double epsrel, int maxit, double min, | ||
899 | double xl, double xu, RMAT(sol)) { | ||
900 | REQUIRES(solr == maxit && solc == 4,BAD_SIZE); | ||
901 | DEBUGMSG("minimize_only_f"); | ||
902 | gsl_function my_func; | ||
903 | my_func.function = only_f_aux_root; | ||
904 | my_func.params = f; | ||
905 | size_t iter = 0; | ||
906 | int status; | ||
907 | const gsl_min_fminimizer_type *T; | ||
908 | gsl_min_fminimizer *s; | ||
909 | // Starting point | ||
910 | switch(method) { | ||
911 | case 0 : {T = gsl_min_fminimizer_goldensection; break; } | ||
912 | case 1 : {T = gsl_min_fminimizer_brent; break; } | ||
913 | case 2 : {T = gsl_min_fminimizer_quad_golden; break; } | ||
914 | default: ERROR(BAD_CODE); | ||
915 | } | ||
916 | s = gsl_min_fminimizer_alloc (T); | ||
917 | gsl_min_fminimizer_set (s, &my_func, min, xl, xu); | ||
918 | do { | ||
919 | double current_min, current_lo, current_hi; | ||
920 | status = gsl_min_fminimizer_iterate (s); | ||
921 | current_min = gsl_min_fminimizer_x_minimum (s); | ||
922 | current_lo = gsl_min_fminimizer_x_lower (s); | ||
923 | current_hi = gsl_min_fminimizer_x_upper (s); | ||
924 | solp[iter*solc] = iter + 1; | ||
925 | solp[iter*solc+1] = current_min; | ||
926 | solp[iter*solc+2] = current_lo; | ||
927 | solp[iter*solc+3] = current_hi; | ||
928 | iter++; | ||
929 | if (status) /* check if solver is stuck */ | ||
930 | break; | ||
931 | |||
932 | status = | ||
933 | gsl_min_test_interval (current_lo, current_hi, 0, epsrel); | ||
934 | } | ||
935 | while (status == GSL_CONTINUE && iter < maxit); | ||
936 | int i; | ||
937 | for (i=iter; i<solr; i++) { | ||
938 | solp[i*solc+0] = iter; | ||
939 | solp[i*solc+1]=0.; | ||
940 | solp[i*solc+2]=0.; | ||
941 | solp[i*solc+3]=0.; | ||
942 | } | ||
943 | gsl_min_fminimizer_free(s); | ||
944 | OK | ||
945 | } | ||
946 | |||
947 | |||
948 | |||
895 | // this version returns info about intermediate steps | 949 | // this version returns info about intermediate steps |
896 | int minimize(int method, double f(int, double*), double tolsize, int maxit, | 950 | int minimize(int method, double f(int, double*), double tolsize, int maxit, |
897 | KRVEC(xi), KRVEC(sz), RMAT(sol)) { | 951 | KRVEC(xi), KRVEC(sz), RMAT(sol)) { |