/******************************************************************************* *
*
* This file is part of the General Hidden Markov Model Library,
* GHMM version __VERSION__, see http://ghmm.org
*
* Filename: ghmm/ghmm/root_finder.c
* Authors: Achim Gaedke
*
* Copyright (C) 1998-2004 Alexander Schliep
* Copyright (C) 1998-2001 ZAIK/ZPR, Universitaet zu Koeln
* Copyright (C) 2002-2004 Max-Planck-Institut fuer Molekulare Genetik,
* Berlin
*
* Contact: schliep@ghmm.org
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Library General Public
* License as published by the Free Software Foundation; either
* version 2 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Library General Public License for more details.
*
* You should have received a copy of the GNU Library General Public
* License along with this library; if not, write to the Free
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*
* This file is version $Revision: 2267 $
* from $Date: 2009-04-24 11:01:58 -0400 (Fri, 24 Apr 2009) $
* last change by $Author: grunau $.
*
*******************************************************************************/
#ifdef HAVE_CONFIG_H
# include "../config.h"
#endif /* */
#ifndef DO_WITH_GSL
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "ghmm_internals.h"
double ghmm_zbrent_AB (double (*func) (double, double, double, double),
double x1, double x2, double tol, double A, double B,
double eps)
{
fprintf (stderr, "Function ghmm_zbrent_AB() not implemented!\n");
exit (1);
/*
double a, b, c;
double fa, fb, fc;
a = min(x1, x2);
fa = (*func)(a);
c = max(x1, x2);
fc = (*func)(c);
b = (c - a)/2.0;
fb = (*func)(b);
while (fabs(c - a) > tol + (tol * min(fabs(a), fabs(c))))
{
r = fb/fc;
s = fb/fa;
t = fa/fc;
p = s * (t * (r - t) * (c - b) - (1.0 - r) * (b - a));
q = (t - 1.0) * (r - 1.0) * (s - 1.0);
x = b + p/q;
if (x > a && x < c)
{
/ * Accept interpolating point * /
if (x < b)
{
c = b;
fc = fb;
b = x;
fb = (*func)(b);
}
else if (x > b)
{
a = b;
fa = fb;
b = x;
fb = (*func)(b);
}
}
else
{
/ * Use bisection * /
}
}
*/
}
#else /* */
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
/* struct for function pointer and parameters except 1st one */
struct parameter_wrapper {
double (*func) (double, double, double, double);
double x2;
double x3;
double x4;
} parameter;
/* calls the given function during gsl root solving iteration
the first parameter variates, all other are kept constant
*/
double function_wrapper (double x, void *p)
{
struct parameter_wrapper *param = (struct parameter_wrapper *) p;
return param->func (x, param->x2, param->x3, param->x4);
}
/*
this interface is used in sreestimate.c
*/
double ghmm_zbrent_AB (double (*func) (double, double, double, double), double x1,
double x2, double tol, double A, double B, double eps)
{
/* initialisation of wrapper structure */
struct parameter_wrapper param;
gsl_function f;
#ifdef HAVE_GSL_INTERVAL /* gsl_interval vanished with version 0.9 */
gsl_interval x;
#endif /* */
gsl_root_fsolver * s;
double tolerance;
int success = 0;
double result = 0;
param.func = func;
param.x2 = A;
param.x3 = B;
param.x4 = eps;
f.function = &function_wrapper;
f.params = (void *) ¶m;
tolerance = tol;
/* initialisation */
#ifdef HAVE_GSL_INTERVAL
# ifdef GSL_ROOT_FSLOVER_ALLOC_WITH_ONE_ARG
x.lower = x1;
x.upper = x2;
s = gsl_root_fsolver_alloc (gsl_root_fsolver_brent);
gsl_root_fsolver_set (s, &f, x);
# else
s = gsl_root_fsolver_alloc (gsl_root_fsolver_brent, &f, x);
# endif
#else /* gsl_interval vanished with version 0.9 */
s = gsl_root_fsolver_alloc (gsl_root_fsolver_brent);
gsl_root_fsolver_set (s, &f, x1, x2);
#endif /* */
/* iteration */
do {
success = gsl_root_fsolver_iterate (s);
if (success == GSL_SUCCESS)
{
#ifdef HAVE_GSL_INTERVAL
gsl_interval new_x;
new_x = gsl_root_fsolver_interval (s);
success = gsl_root_test_interval (new_x, tolerance, tolerance);
#else /* gsl_interval vanished with version 0.9 */
double x_up;
double x_low;
(void) gsl_root_fsolver_iterate (s);
x_up = gsl_root_fsolver_x_upper (s);
x_low = gsl_root_fsolver_x_lower (s);
success = gsl_root_test_interval (x_low, x_up, tolerance, tolerance);
#endif /* */
}
} while (success == GSL_CONTINUE);
/* result */
if (success != GSL_SUCCESS)
{
gsl_error ("solver failed", __FILE__, __LINE__, success);
}
else
{
result = gsl_root_fsolver_root (s);
}
/* destruction */
gsl_root_fsolver_free (s);
return result;
}
#endif /* */