summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--lib/Data/Packed/Development.hs1
-rw-r--r--lib/Numeric/GSL/Minimization.hs34
-rw-r--r--lib/Numeric/GSL/gsl-aux.c54
3 files changed, 89 insertions, 0 deletions
diff --git a/lib/Data/Packed/Development.hs b/lib/Data/Packed/Development.hs
index 9b8dacf..e792783 100644
--- a/lib/Data/Packed/Development.hs
+++ b/lib/Data/Packed/Development.hs
@@ -21,6 +21,7 @@ module Data.Packed.Development (
21 app1, app2, app3, app4, 21 app1, app2, app3, app4,
22 app5, app6, app7, app8, app9, app10, 22 app5, app6, app7, app8, app9, app10,
23 MatrixOrder(..), orderOf, cmat, fmat, 23 MatrixOrder(..), orderOf, cmat, fmat,
24 matrixFromVector,
24 unsafeFromForeignPtr, 25 unsafeFromForeignPtr,
25 unsafeToForeignPtr, 26 unsafeToForeignPtr,
26 check, (//), 27 check, (//),
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 756edf1..24d82c4 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>
@@ -905,6 +906,59 @@ double only_f_aux_min(const gsl_vector*x, void *pars) {
905 return res; 906 return res;
906} 907}
907 908
909double only_f_aux_root(double x, void *pars);
910int uniMinimize(int method, double f(double),
911 double epsrel, int maxit, double min,
912 double xl, double xu, RMAT(sol)) {
913 REQUIRES(solr == maxit && solc == 4,BAD_SIZE);
914 DEBUGMSG("minimize_only_f");
915 gsl_function my_func;
916 my_func.function = only_f_aux_root;
917 my_func.params = f;
918 size_t iter = 0;
919 int status;
920 const gsl_min_fminimizer_type *T;
921 gsl_min_fminimizer *s;
922 // Starting point
923 switch(method) {
924 case 0 : {T = gsl_min_fminimizer_goldensection; break; }
925 case 1 : {T = gsl_min_fminimizer_brent; break; }
926 case 2 : {T = gsl_min_fminimizer_quad_golden; break; }
927 default: ERROR(BAD_CODE);
928 }
929 s = gsl_min_fminimizer_alloc (T);
930 gsl_min_fminimizer_set (s, &my_func, min, xl, xu);
931 do {
932 double current_min, current_lo, current_hi;
933 status = gsl_min_fminimizer_iterate (s);
934 current_min = gsl_min_fminimizer_x_minimum (s);
935 current_lo = gsl_min_fminimizer_x_lower (s);
936 current_hi = gsl_min_fminimizer_x_upper (s);
937 solp[iter*solc] = iter + 1;
938 solp[iter*solc+1] = current_min;
939 solp[iter*solc+2] = current_lo;
940 solp[iter*solc+3] = current_hi;
941 iter++;
942 if (status) /* check if solver is stuck */
943 break;
944
945 status =
946 gsl_min_test_interval (current_lo, current_hi, 0, epsrel);
947 }
948 while (status == GSL_CONTINUE && iter < maxit);
949 int i;
950 for (i=iter; i<solr; i++) {
951 solp[i*solc+0] = iter;
952 solp[i*solc+1]=0.;
953 solp[i*solc+2]=0.;
954 solp[i*solc+3]=0.;
955 }
956 gsl_min_fminimizer_free(s);
957 OK
958}
959
960
961
908// this version returns info about intermediate steps 962// this version returns info about intermediate steps
909int minimize(int method, double f(int, double*), double tolsize, int maxit, 963int minimize(int method, double f(int, double*), double tolsize, int maxit,
910 KRVEC(xi), KRVEC(sz), RMAT(sol)) { 964 KRVEC(xi), KRVEC(sz), RMAT(sol)) {