root_finder.c

来自「General Hidden Markov Model Library 一个通用」· C语言 代码 · 共 217 行

/*******************************************************************************  *
*
*       This file is part of the General Hidden Markov Model Library,
*       GHMM version 0.8_beta1, 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: [e-mail]
*
*       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: 1439 $
*                       from $Date: 2005-10-12 16:35:51 +0200 (Wed, 12 Oct 2005) $
*             last change by $Author: grunau $.
*
*******************************************************************************/
   
  
#ifdef WIN32
#  include "win_config.h"
#endif  /*  */
  
#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 *) &param;
  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  /*  */