COPASI API  4.16.103
Classes | Functions
FminBrent.h File Reference
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes

class  FDescent
 
class  FDescentTemplate< CType >
 

Functions

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

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