File: fmin.i

package info (click to toggle)
yorick-mira 0.9.9+dfsg1-2
  • links: PTS
  • area: main
  • in suites: squeeze
  • size: 1,280 kB
  • ctags: 3
  • sloc: makefile: 90
file content (276 lines) | stat: -rw-r--r-- 9,591 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
/*
 * fmin.i -
 *
 *	Minimization of an univariate function for Yorick.  The method is
 *	based on original Brent's method modified to allow for different
 *	kind of bounds (both, left, right or none).
 *
 *-----------------------------------------------------------------------------
 *
 *      Copyright (C) 2001, Eric Thi´┐Żbaut <thiebaut@obs.univ-lyon1.fr>
 *
 *	This file is free software; you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License version 2 as
 *	published by the Free Software Foundation.
 *
 *	This file 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.
 *
 *-----------------------------------------------------------------------------
 *
 *	$Id: fmin.i,v 1.2 2008/07/12 06:47:24 eric Exp $
 *	$Log: fmin.i,v $
 *	Revision 1.2  2008/07/12 06:47:24  eric
 *	 - Added final comment for setting local variables of Emacs.
 *
 *	Revision 1.1  2005/06/24 18:14:35  eric
 *	Initial revision
 *
 *-----------------------------------------------------------------------------
 */

func fmin(f, a, b, lim, tol=, all=, eps=)
/* DOCUMENT fmin(f, a, b)
       -or- fmin(f, a, b, lim)
     Get the location  of the minimum of univariate function  F(X).  F is a
     Yorick  function,  A and  B  are bounds  or  starting  values for  the
     variable X and optional LIM specifies the kind of limits for X:
       If LIM=0 or nil, there is no bounds for X: F is first evaluated at A
         and B, then the interval of  search is enlarged until a minimum is
         bracketed.
       If LIM=1,  the interval is bounded by  A: F is first  evaluated at B
         and  the interval is  enlarged (away  from A)  until a  minimum is
         bracketed  --  i.e.  the location of  the  minimum X is such that:
         A < X, if A < B; or A > X, if A > B.
       If LIM=2, the  interval is bounded by B (same as  with LIM=1 but the
         role of A and B exchanged).
       If LIM=3, the minimum is searched  in the interval (A,B) -- i.e. the
         location of the minimum X is such that: min(A,B) < X < max(A,B).

     Keyword  TOL can be  used to  specify the  relative precision  for the
     solution.  The default value is TOL=sqrt(EPS) (see below).

     If keyword ALL  is true (non-nil and non-zero)  the returned value is:
     [X, FX, XLO, XHI] where X is the approximated location of the minimum,
     FX is F(X) and XLO and XHI are the lower and upper bounds for the true
     minimum; otherwise, only X is returned.

     Keyword EPS  can be  used to specify  the machine  relative precision.
     Default value is EPS~2.22e-16,  which corresponds to IEEE standard for
     double precision.

   NOTES:
     (1) The  minimization routine  should never evaluates  F at  the given
         bounds if any.
     (2) If the function F(X) is  not unimodal, only a local minimum can be
         found.

   REFERENCES:
     The method is  based on original Brent's method  modified to allow for
     different kind of bounds (both, left, right or none).

   SEE ALSO: */
{
  /* Make sure A and B are double precision values. */
  a += 0.0;
  b += 0.0;

  /* EPS is approximately the square root of the relative machine
     precision. */
  if (is_void(eps)) eps = 2.2204460492503131e-16; /* assume IEEE double */
  tol1 = eps + 1.0;
  eps = sqrt(eps);
  tol3 = (is_void(tol) ? eps : tol)/3.0;
  /* TOL not used below */

  /* C = (3 - sqrt(5))/2 is the squared inverse of the golden ratio */
  c = 0.3819660112501051517954131656343618822796908201942371378645513772947395;

  /* S = (1 + sqrt(5))/2 = 2 - C is a constant used to increase the width
     of the interval with a golden ratio. */
  s = 1.6180339887498948482045868343656381177203091798057628621354486227052605;

  /* Original Brent's method assumes that the minimum is in (A,B) with
   * A<=B and keeps track of the following variables:
   *   X, FX = least function value found so far
   *   W, FW = previous value of X, FX
   *   V, FV = previous value of W, FW
   *   U, FU = last function evaluation
   * If the interval to consider is not bounded or only left/right bounded,
   * the idea is to find a suitable interval (A,B) where at least one
   * minimum must exists (if the function is continue) and start Brent's
   * algorithm with correct values for X, FX, ... (in order to save some
   * function evaluations).
   */
  if (! lim) {
    /* The interval of search is unlimited, we start with A, B and then
       search for a bracket. */
    x = a;
    fx = f(x);
    w = b;
    fw = f(w);

    /* Make sure X is the best location found so far, we therefore exchange
       W and X if FW <= FX (we exchange the two point in case of equality
       to alternatively search on the other side where the function is,
       numerically, flat) */
    if (fw <= fx) {
      tmp=w; w=x; x=tmp;
      tmp=fw; fw=fx; fx=tmp;
    }

    /* Loop until a bracket is found. Possible improvements: (1) use parabolic
       extrapolation to allow for bigger jumps, (2) keep track of one more
       point in order to be able to slightly reduce the size of the interval
       once a bracket has been found. */
    for (;;) {
      /* Take a golden step in the descent direction. */
      v = x + s*(x - w);
      fv = f(v);

      if (fw > fx) {
        if (fv > fx) {
          /* Bracket found: the minimum is in (V,W).  Set variables for
             Brent's method: set bounds such that A is smaller than B. */
          if (v < w) {
            a = v;
            b = w;
          } else {
            a = w;
            b = v;
          }
          break; /* branch to Brent's method */
        } else {
          /* Continue with golden search with X and V. */
          w=x; fw=fx;
          x=v; fx=fv;
        }
      } else {
        /* We are in a, numerically, flat region (FW=FX). Enlarge interval
           (V,W) by a factor 1+S ~ 2.62 */
        if (fv >= fx) {
          x=w; fx=fw;
          w=v; fw=fv;
        } else {
          w=x; fw=fx;
          x=v; fx=fv;
        }
      }
    }
  } else if (lim == 1 || lim == 2) {
    /* Interval is bounded by A or B.  Possibly exchange A and B, so that
       A is the bound and search until a bound for the other side is found. */
    if (lim == 2) {tmp=a; a=b; b=tmp;}
    w = x = b;
    fw = fx = f(x);
    for (;;) {
      v = x + s*(x - a);
      fv = f(v);
      if (fv > fx) {
        /* We have found a bound for the other side.  Set search interval
           to be (A,B) := (A,V) or (V,A) such that A is smaller than B. */
        if (v > a) {
          b = v;
        } else {
          b = a;
          a = v;
        }
        break; /* branch to Brent's method */
      }
      w=x; fw=fx;
      x=v; fx=fv;
      if (fw > fx) a = w;
    }
  } else if (lim == 3) {
    /* The minimum is to be found in (A,B) -- this is original Brent's
       method.  Make sure that A is smaller than B and set X,FX ... */
    if (a > b) {tmp=a; a=b; b=tmp;}
    v = w = x = a + c*(b - a);
    fv = fw = fx = f(x);
  } else {
    error, "bad value for keyword LIM";
  }

  /*** Brent's method. ***/

  /* Set E and D (note: the golden step instead of the parabolic step is
     taken if abs(E) is too small). */
  e = x - v;
  d = x - w;

  /* main loop starts here */
  for (;;) {
    xm = (a + b)*0.5;
    tol1 = eps*abs(x) + tol3;
    tol2 = tol1 + tol1;

    /* check stopping criterion */
    if (abs(x - xm) <= tol2 - (b - a)*0.5) {
      if (all) return [x, fx, a, b];
      return x;
    }
    if (abs(e) > tol1) {
      /* fit parabola */
      q = (x - v)*(fx - fw);
      r = (x - w)*(fx - fv);
      if (q <= r) {
	p = (x - v)*q - (x - w)*r;
	q = (r - q)*2.0;
      } else {
	p = (x - w)*r - (x - v)*q;
	q = (q - r)*2.0;
      }
      if (abs(p) < abs(0.5*q*e) && p > q*(a - x) && p < q*(b - x)) {
	/* use a parabolic-interpolation step */
        e = d;
        u = x + (d = p/q);

	/* F must not be evaluated too close to A or B */
	if (u - a < tol2 || b - u < tol2) {
	  d = (x < xm ? tol1 : -tol1);
	}
      } else {
	/* use a golden-section step */
	e = (x >= xm ? a : b) - x;
	d = c*e;
      }
    } else {
      /* use a golden-section step */
      e = (x >= xm ? a : b) - x;
      d = c*e;
    }

    /* F must not be evaluated too close to X */
    u = (abs(d) >= tol1 ? x + d : (d > 0.0 ? x + tol1 : x - tol1));
    fu = f(u);

    /* update A, B, V, W, and X */
    if (fx <= fu) {
      if (u >= x) b = u;
      else        a = u;
    }
    if (fu <= fx) {
      if (u >= x) a = x;
      else        b = x;
      v = w; fv = fw;
      w = x; fw = fx;
      x = u; fx = fu;
    } else if (fu <= fw || w == x) {
      v = w; fv = fw;
      w = u; fw = fu;
    } else if (fu <= fv || v == x || v == w) {
      v = u; fv = fu;
    }
  } /* end of main loop */
}

/*---------------------------------------------------------------------------*
 * Local Variables:                                                          *
 * mode: Yorick                                                              *
 * tab-width: 8                                                              *
 * fill-column: 75                                                           *
 * c-basic-offset: 2                                                         *
 * coding: latin-1                                                           *
 * End:                                                                      *
 *---------------------------------------------------------------------------*/