Commit 63a1660b authored by zhanxw's avatar zhanxw
Browse files

Improve GSLMinimizer robustness

parent 48dc8723
......@@ -9,18 +9,28 @@ Minimizer::Minimizer()
T = gsl_min_fminimizer_brent;
s = gsl_min_fminimizer_alloc(T);
}
Minimizer::~Minimizer() { gsl_min_fminimizer_free(s); }
int Minimizer::minimize(gsl_function F, double startValue, double lowerBound,
double upperBound) {
gsl_min_fminimizer_set(s, &F, startValue, lowerBound, upperBound);
gsl_error_handler_t* oldHandler = gsl_set_error_handler_off();
int iter = 0;
int status;
double a, b;
status = gsl_min_fminimizer_set(s, &F, startValue, lowerBound, upperBound);
if (status != GSL_SUCCESS) {
#ifndef NDEBUG
fprintf(stderr, "Minimizer failed due to: %s\n", gsl_strerror(status));
#endif
goto errorLabel;
}
do {
iter++;
status = gsl_min_fminimizer_iterate(s);
if (status == GSL_EBADFUNC || status == GSL_FAILURE) {
return -1;
goto errorLabel;
}
finalX = gsl_min_fminimizer_x_minimum(s);
......@@ -35,7 +45,7 @@ int Minimizer::minimize(gsl_function F, double startValue, double lowerBound,
if (status == GSL_SUCCESS) {
// printf ("Converged:\n");
finalY = gsl_min_fminimizer_f_minimum(s);
return 0;
goto successLabel;
}
// printf("%5d\ta=%g\tb=%g\tfinalX=%g\n", iter, a, b, finalX);
// printf ("%5d .7f, [.7f] "
......@@ -44,5 +54,10 @@ int Minimizer::minimize(gsl_function F, double startValue, double lowerBound,
// m, m - m_expected, b - a);
} while (status == GSL_CONTINUE && iter < this->maxIter);
successLabel:
gsl_set_error_handler(oldHandler);
return 0;
errorLabel:
gsl_set_error_handler(oldHandler);
return -1;
}
......@@ -28,6 +28,7 @@ class Minimizer {
// GSL stuffs
const gsl_min_fminimizer_type *T;
gsl_min_fminimizer *s;
int status;
};
#endif /* _GSLMINIMIZER_H_ */
......@@ -15,6 +15,7 @@ EXE = testGSLIntegration \
testFastMultipleTraitScoreTest \
testFastMultipleTraitScoreTest.NA \
testBoltLMM \
testGSLMinimizer \
all: $(EXE)
debug: $(EXE)
......
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include "regression/GSLMinimizer.h"
#include "third/gsl/include/gsl/gsl_errno.h"
static double func1(double x, void* param) { return x * x; }
static double func2(double x, void* param) {
// this example is taken from:
// https://stackoverflow.com/questions/18807512/is-this-a-bug-in-gsls-minimum-finding-routine
// direct optimization will make it fail
double beta = x;
return -(
((6160558822864 * (exp(4 * beta)) + 523830424923 * (exp(3 * beta)) +
1415357447750 * (exp(5 * beta)) + 7106224104 * (exp(6 * beta))) /
(385034926429 * (exp(4 * beta)) + 58203380547 * (exp(3 * beta)) +
56614297910 * (exp(5 * beta)) + 197395114 * (exp(6 * beta)))) -
((1540139705716 * (exp(4 * beta)) + 174610141641 * (exp(3 * beta)) +
283071489550 * (exp(5 * beta)) + 1184370684 * (exp(6 * beta))) *
(1540139705716 * (exp(4 * beta)) + 174610141641 * (exp(3 * beta)) +
283071489550 * (exp(5 * beta)) + 1184370684 * (exp(6 * beta))) /
pow((385034926429 * (exp(4 * beta)) + 58203380547 * (exp(3 * beta)) +
56614297910 * (exp(5 * beta)) + 197395114 * (exp(6 * beta))),
2)));
}
int main() {
{
Minimizer m;
gsl_function F;
F.function = func1;
F.params = NULL;
double start = 0.5;
double lb = -5;
double ub = 5;
if (m.minimize(F, start, lb, ub)) {
fprintf(stderr, "Minimizer failed");
assert(false);
} else {
fprintf(stderr, "Minimizer succeed, x = %g\n", m.getX());
}
}
{
Minimizer m;
gsl_function F;
F.function = func2;
F.params = NULL;
double start = 0.0;
double lb = -6;
double ub = 6;
// gsl_set_error_handler_off();
if (m.minimize(F, start, lb, ub)) {
fprintf(stderr, "Minimizer failed\n");
} else {
fprintf(stderr, "Minimizer succeed, x = %g\n", m.getX());
assert(false);
}
}
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment