COPASI API  4.16.103
FminBrent.cpp File Reference
`#include <cmath>`
`#include <limits>`
`#include "copasi.h"`
`#include "FminBrent.h"`
Include dependency graph for FminBrent.cpp:

Go to the source code of this file.

## Macros

#define SQRT_EPSILON   sqrt(std::numeric_limits< C_FLOAT64 >::epsilon())

## Functions

int FminBrent (double a, double b, FDescent *pF, double *min, double *fmin, double tol, int maxiter)

## Macro Definition Documentation

 #define SQRT_EPSILON   sqrt(std::numeric_limits< C_FLOAT64 >::epsilon())

Definition at line 29 of file FminBrent.cpp.

Referenced by FminBrent().

## Function Documentation

 int FminBrent ( double a, double b, FDescent * pF, double * min, double * fmin, double tol, int maxiter )

Definition at line 30 of file FminBrent.cpp.

References SQRT_EPSILON.

Referenced by COptMethodSteepestDescent::optimise().

37 {
38  double x, v, w; /* Abscissae, descr. see above */
39  double fx; /* f(x) */
40  double fv; /* f(v) */
41  double fw; /* f(w) */
42  const double r = (3. - sqrt(5.0)) / 2; /* Gold section ratio */
43  int iter; /* Iteration counter */
44
45  if (tol <= 0) return 1; /* check input values */
46
47  if (b <= a) return 2;
48
49  v = a + r * (b - a); fv = (*pF)(v); /* First step - always gold section*/
50  x = v; w = v;
51  fx = fv; fw = fv;
52
53  for (iter = 0; iter < maxiter; iter++) /* Main iteration loop */
54  {
55  double range = b - a; /* Range over which the minimum is */
56  double middle_range = (a + b) / 2; /* seeked */
57  double tol_act = /* Actual tolerance */
58  SQRT_EPSILON * fabs(x) + tol / 3;
59  double new_step; /* Step at this iteration */
60
61  if (fabs(x - middle_range) + range / 2 <= 2*tol_act)
62  {
63  *min = x; *fmin = fx; /* Store the solution */
64  return 0; /* Acceptable approx. is found */
65  }
66
67  /* Obtain the gold section step */
68  new_step = r * (x < middle_range ? b - x : a - x);
69
70  /* Decide if the interpolation */
71
72  /* can be tried */
73  if (fabs(x - w) >= tol_act) /* If x and w are distinct */
74  {/* interpolatiom may be tried */
75  double p; /* Interpolation step is calculated */
76  double q; /* as p/q; division operation */
77  double t; /* is delayed until last moment */
78
79  t = (x - w) * (fx - fv);
80  q = (x - v) * (fx - fw);
81  p = (x - v) * q - (x - w) * t;
82  q = 2 * (q - t);
83
84  if (q > 0.0) /* q was calculated with the */
85  p = -p; /* opposite sign; make q positive */
86  else /* and assign possible minus to */
87  q = -q; /* p */
88
89  if (fabs(p) < fabs(new_step*q) && /* If x+p/q falls in [a,b] */
90  p > q*(a - x + 2*tol_act) && /* not too close to a and */
91  p < q*(b - x - 2*tol_act)) /* b, and isn't too large */
92  new_step = p / q; /* it is accepted */
93
94  /* If p/q is too large then the */
95  /* gold section procedure can */
96  /* reduce [a,b] range to more */
97  /* extent */
98  }
99
100  if (fabs(new_step) < tol_act) /* Adjust the step to be not less */
101  {
102  if (new_step > (double)0) /* than tolerance */
103  new_step = tol_act;
104  else
105  new_step = -tol_act;
106  }
107
108  /* Obtain the next approximation to */
109  {/* min & reduce the enveloping range*/
110  double t = x + new_step; /* Tentative point for the min */
111  double ft = (*pF)(t);
112
113  if (ft <= fx)
114  {/* t is a better approximation */
115  if (t < x) /* Reduce the range so that */
116  b = x; /* t would fall within it */
117  else
118  a = x;
119
120  v = w; w = x; x = t; /* Assign the best approx to x */
121  fv = fw; fw = fx; fx = ft;
122  }
123  else /* x remains the better approx */
124  {
125  if (t < x) /* Reduce the range enclosing x */
126  a = t;
127  else
128  b = t;
129
130  if (ft <= fw || w == x)
131  {
132  v = w; w = t;
133  fv = fw; fw = ft;
134  }
135  else if (ft <= fv || v == x || v == w)
136  {
137  v = t;
138  fv = ft;
139  }
140  }
141  } /* ----- end-of-block ----- */
142  } /* ===== End of loop ===== */
143
144  *min = x; *fmin = fx; /* Store the best value */
145  return 3; /* Too many iterations */
146 }
#define SQRT_EPSILON
Definition: FminBrent.cpp:29
#define min(a, b)
Definition: f2c.h:175