summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--lib/Numeric/GSL/Minimization.hs34
-rw-r--r--lib/Numeric/GSL/gsl-aux.c54
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
53module Numeric.GSL.Minimization ( 53module 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
85data UniMinimizeMethod = GoldenSection
86 | BrentMini
87 | QuadGolden
88 deriving (Enum, Eq, Show, Bounded)
89
90-- | Onedimensional minimization.
91
92uniMinimize :: 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
101uniMinimize method epsrel maxit fun xmin xl xu = uniMinimizeGen (fi (fromEnum method)) fun xmin xl xu epsrel maxit
102
103uniMinimizeGen 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
115foreign import ccall safe "uniMinimize"
116 c_uniMinize:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> Double -> TM
117
84data MinimizeMethod = NMSimplex 118data 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
896double only_f_aux_root(double x, void *pars);
897int 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
896int minimize(int method, double f(int, double*), double tolsize, int maxit, 950int 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)) {