Integrated Search using Parabolic, Bracket, and Downhill
Language: English
Programming Language: C
Published by: n_nelson [not registered]
Last Update: 5/15/2006
Views: 148
Description
This routine provides an optimization (on y, minima or maxima) search of a function with one independent variable, y=f(x). The three search methods of parabolic, bracket, and downhill are integrated in that prioirity sequence. The method performs only the optimization function external to the function being searched, keeps separate the search state variables, and hence may be used for simultaneous parallel or searches within searches. Several optional search parameters are provided.
Code
1 /* PSRCH.C ver: 0.1 - Parabolic search routine for a single dimension.
2
3 This code may be copied and used without restriction of any kind.
4 No warranties are given or implied.
5
6 Author: Neil Nelson <n_nelson@pacbell.net>, February, 21, 2001.
7
8 This is a more robust version of my previous psearch.c.
9
10 * A search proceeds first via a parabolic estimation, then if for some reason the parabolic
11 estimation is not acceptable (out of range or inverted), the search uses a bracket search
12 when high and low bounds have been identified via any previous estimation, and then if
13 both brackets are not present, the search uses a downhill method.
14
15 * The previous psearch.c and those shown in _Numerical Recipes in C_ are passed a function
16 address, and then subsequent activity for the function that is searched is embedded
17 within or below the control of the search routine. I wanted greater direct access to the
18 routine being searched and to be able to have the capability of having a number of
19 searches going on simultaneously either in parallel or searches within searches. This
20 resulted in a routine that performs on each call one cycle of the search process and
21 utilizes on each call the state variables for that particular search. Note that there
22 will commonly be a minimum of four calls to psrch, the first three to set up the
23 parabola and then a fourth call to check, at a minimum, the first estimate.
24
25 * The include file psrch.h for the state variables is appended below.
26
27 *** Modifications
28 * 010225 (1) Implemented max_jmp.
29 (2) Corrected pt return on search done. Prior return did not match the yr return
30 and was not the best of those available.
31
32 Example calling program:
33 psrch_vars *psrch_id; // Pointer to search parameters and work values. All the
34 // state variables for the search are held here so that the
35 // psrch routines may be used for a variety of searches
36 // at the same time.
37 double pt,y=0.0; // Search start pt and a variable to pass the value of the
38 // function searched. y is to minimized or maximized
39
40 psrch_vars *psrch_open(); // Set up state variables.
41 int psrch_getpt(psrch_vars *psrch_id, double *pt, double *y); // Optimization routine.
42 void psrch_close(psrch_vars *psrch_id); // Free state variables.
43 double function_searched(double pt); // A function that takes pt as an independent
44 // variable and returns y or however you wish to associate pt with y.
45
46 psrch_id=psrch_open(); // Obtains pointer to state variables.
47 if (!psrch_id) return 0;
48
49 psrch_id->srch_typ=(MAXIMA); // Search type: 1=MAXIMA, 0=MINIMA
50 psrch_id->srch_lmt=0.0001; // Search ends when the difference between the new y and
51 // the current best y of the search is less than srch_lmt.
52
53 pt=??; // Set the value of starting point from which the search will begin.
54
55 psrch_id->dflt=1.0; // Default increment for pt. Used to get direction of search.
56 psrch_id->units=0; // Is this a unit/integer dimension? (Unit searches not well tested.)
57
58 // ** Optional parameters that are defaulted but may be set by the user.
59 // psrch_id->max_jmp // (double) The largest search increment that may be used
60 // for a change in pt.
61 // Defaults to around 10 times srch_lmt.
62 // psrch_id->imprv_fctr_lmt // (float) If imprv_fctr is greater than imprv_fctr_lmt,
63 // then done. Defaults to 3.0;
64 // psrch_id->prtl_imprv; // (float) Amount added to imprv_fctr when less than best
65 // improvement obtained. When this value goes over
66 // imprv_fctr_lmt search is done. Default is 0.34.
67 // psrch_id->brkt_lmt; // (int) Max allowed number of converging bracket pairs
68 // without better y.
69
70 // ** Options not yet implemented
71 // int lmt_cnfms; // NA When gain limit reached, how many tries to get greater
72 // than srch_lmt.
73 // double max_hi; // NA Largest value allowed.
74 // double max_lo; // NA Smallest value allowed.
75
76 for(i=0;i<20;i++) { // Loop for maximum number of search passes. Should be greater than 3.
77 if (psrch_getpt(psrch_id, &pt, &y)) break; // Search done.
78 y=function_searched(pt); // else try next pt.
79 }
80
81 */
82
83 #include <stdio.h>
84 #include <stdlib.h>
85 #include <math.h>
86 #include "psrch.h"
87
88 void psrch_close(psrch_vars *psrch_id)
89 {
90 if (psrch_id) {
91 free(psrch_id);
92 }
93 return;
94 }
95
96 psrch_vars *psrch_open() {
97 psrch_vars *psrch_id;
98
99 psrch_id = (psrch_vars*) malloc(sizeof(psrch_vars));
100 if (!psrch_id) {
101 printf("%s\n", "psrch_open - psrch_vars malloc error");
102 return psrch_id;
103 }
104
105 psrch_id->dflt=0.0;
106 psrch_id->units=0;
107 psrch_id->strtpt=0.0;
108 psrch_id->ptsnt=0.0;
109 psrch_id->max_hi=0.0;
110 psrch_id->max_lo=0.0;
111
112 psrch_id->srch_typ=0; // Minima search.
113 psrch_id->srch_fnc=0; // Parabolic search.
114 psrch_id->max_jmp=0.0;
115 psrch_id->srch_lmt=0.0;
116 psrch_id->lmt_cnfms=3;
117 psrch_id->imprv_fctr=0.0;
118 psrch_id->imprv_fctr_lmt=3.0;
119 psrch_id->prtl_imprv=0.34;
120 psrch_id->ygain=0; // Not yet any gain in y.
121 psrch_id->srch_hix=psrch_id->srch_lox=0.0;
122 psrch_id->brkset=psrch_id->brksde=0; // Bracket flag set to neither side.
123 psrch_id->strty=0.0;
124 psrch_id->xlngsnt=0.0;
125 psrch_id->call_code=-1;
126 psrch_id->brkt_dvsr=1.0;
127 psrch_id->brkt_cnt=0;
128 psrch_id->brkt_lmt=2;
129 psrch_id->ysrt[0]=0;
130
131 return psrch_id;
132 }
133
134 void psrch_srt(psrch_vars *psrch_id) {
135 /*
136 * This routine obtains an ascending (by y_seq index) order
137 * of the current values of the dependent (y) variable.
138 */
139 int *ysrt;
140 double *y;
141
142 y=psrch_id->y;
143 ysrt=psrch_id->ysrt;
144
145 if (y[0] > y[1]) {
146 ysrt[2] = 0;
147 ysrt[1] = 1;
148 }
149 else {
150 ysrt[2] = 1;
151 ysrt[1] = 0;
152 }
153
154 if (y[2] > y[ysrt[2]]) {
155 ysrt[0] = ysrt[1];
156 ysrt[1] = ysrt[2];
157 ysrt[2] = 2;
158 }
159 else if (y[2] > y[ysrt[1]]) {
160 ysrt[0] = ysrt[1];
161 ysrt[1] = 2;
162 }
163 else ysrt[0] = 2;
164
165 return;
166 }
167
168 int psrch_brckt(psrch_vars *psrch_id) { // Bracket method
169
170 if ((psrch_id->brkt_cnt/2) <= psrch_id->brkt_lmt) {
171 psrch_id->brkt_cnt += 1;
172 if (!psrch_id->brksde) psrch_id->brkt_dvsr /= BRACKET_FACTOR;
173
174 psrch_id->xlngsnt = psrch_id->xlng[psrch_id->ysrt[0]]; // start at best pt
175 // Use opposite side of last xlngsnt bracket first.
176 if (psrch_id->brksde != 1
177 && psrch_id->xlngsnt > psrch_id->xlng[psrch_id->ysrt[0]]) {
178 psrch_id->xlngsnt += -(fabs(psrch_id->xlng[psrch_id->ysrt[0]]
179 - psrch_id->srch_lox)) / psrch_id->brkt_dvsr;
180 psrch_id->brksde=(!psrch_id->brksde ? 1 : 0);
181
182 } else {
183 psrch_id->xlngsnt += (fabs(psrch_id->xlng[psrch_id->ysrt[0]]
184 - psrch_id->srch_hix)) / psrch_id->brkt_dvsr;
185 psrch_id->brksde=(!psrch_id->brksde ? 2 : 0);
186 }
187 return 1;
188 } else return 0;
189 }
190
191 void psrch_cbrckt(psrch_vars *psrch_id, double y, double xlnsnt) {
192
193 if (!psrch_id->brkset) {
194 psrch_id->brkt_cnt=psrch_id->brksde=0;
195 psrch_id->brkt_dvsr=1.0;
196 psrch_id->srch_lox=psrch_id->srch_hix=0.0;
197 }
198
199 // Get brackets for jump boundaries and linear search.
200 if ( (xlnsnt > psrch_id->srch_lox || psrch_id->srch_lox == 0.0)
201 && xlnsnt < psrch_id->xlng[psrch_id->ysrt[0]]
202 && y > psrch_id->y[psrch_id->ysrt[0]]) {
203 psrch_id->srch_lox = xlnsnt;
204
205 if (psrch_id->brkset != 3) {
206 if (psrch_id->brkset == 0) psrch_id->brkset=1;
207 else if (psrch_id->brkset == 2) psrch_id->brkset=3;
208 }
209 return;
210
211 } else if ( (xlnsnt < psrch_id->srch_hix || psrch_id->srch_hix == 0.0)
212 && xlnsnt > psrch_id->xlng[psrch_id->ysrt[0]]
213 && y > psrch_id->y[psrch_id->ysrt[0]]) {
214 psrch_id->srch_hix = xlnsnt;
215
216 if (psrch_id->brkset != 3) {
217 if (psrch_id->brkset == 0) psrch_id->brkset=2;
218 else if (psrch_id->brkset == 1) psrch_id->brkset=3;
219 }
220 return;
221 }
222
223 return;
224 }
225
226 void psrch_dwnhll(psrch_vars *psrch_id) {
227 /* Checks OK.
228 Determines an xlngsnt (pt) that is on the down hill side of the best yet found pt.
229 */
230 int *ysrt;
231 double td1,td2;
232
233 ysrt=psrch_id->ysrt;
234
235 td1 = psrch_id->xlng[ysrt[0]] - psrch_id->xlng[ysrt[2]];
236 td2 = psrch_id->xlng[ysrt[0]] - psrch_id->xlng[ysrt[1]];
237
238 if (fabs((double)td1) < fabs((double)td2))
239 psrch_id->xlngsnt = psrch_id->xlng[ysrt[0]] + td1 / INVERSION_DIVISOR;
240 else psrch_id->xlngsnt = psrch_id->xlng[ysrt[0]] + td2 / INVERSION_DIVISOR;
241
242 td1 = psrch_id->xlngsnt - psrch_id->xlng[ysrt[0]];
243 if (fabs(td1) < TINY)
244 psrch_id->xlngsnt=(td1 < 0
245 ? psrch_id->xlng[ysrt[0]]-TINY : psrch_id->xlng[ysrt[0]]+TINY);
246
247 if (fabs(psrch_id->xlngsnt - psrch_id->xlng[ysrt[0]]) > psrch_id->max_jmp)
248 psrch_id->xlngsnt=((psrch_id->xlngsnt - psrch_id->xlng[ysrt[0]]) > 0
249 ? psrch_id->xlng[ysrt[0]]+psrch_id->max_jmp
250 : psrch_id->xlng[ysrt[0]]-psrch_id->max_jmp);
251 }
252
253 int psrch_parest(psrch_vars *psrch_id) {
254 int *ysrt,prb_ok=1;
255 double yi[3],*xlng,a,b=0.0,*y;
256
257 y=psrch_id->y;
258 xlng=psrch_id->xlng;
259 ysrt=psrch_id->ysrt;
260 xlng=psrch_id->xlng;
261
262 // Alternate parabolic minimum compute from Numerical Recipes.
263 // r=(*bx-*ax)*(*fb-*fc); // Compute u by parabolic extrapolation from a,b,c.
264 // q=(*bx-*cx)*(*fb-*fa); // TINY is used to prevent any possible division by zero.
265 // u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/(2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));
266
267 // Compute intermediate Lagrangian interpolation values.
268 yi[0] = y[0] / ((xlng[0] - xlng[1]) * (xlng[0] - xlng[2]));
269 yi[1] = y[1] / ((xlng[1] - xlng[0]) * (xlng[1] - xlng[2]));
270 yi[2] = y[2] / ((xlng[2] - xlng[0]) * (xlng[2] - xlng[1]));
271
272 // Compute parabolic equation coeficients: ax^2+bx+c.
273 a = yi[0] + yi[1] + yi[2];
274 if (fabs((double)a) < TINY) a=(a > 0.0 ? TINY : -TINY);
275
276 if (a < 0) { // Inverted parabolic curve.
277 prb_ok=0;
278
279 } else {
280
281 if (fabs((double)a) > LARGE) {
282 if (a > 0.0) psrch_id->xlngsnt=(-b > 0.0 ? TINY : -TINY);
283 else psrch_id->xlngsnt=(b > 0.0 ? TINY : -TINY);
284
285 } else {
286 b = - ( yi[0] * (xlng[1] + xlng[2])
287 + yi[1] * (xlng[0] + xlng[2])
288 + yi[2] * (xlng[0] + xlng[1]));
289
290 psrch_id->xlngsnt = - b / (2.0 * a); // maxima or minima of parabolic equation.
291 if (fabs((psrch_id->xlngsnt - xlng[ysrt[0]])/xlng[ysrt[0]])
292 < psrch_id->srch_lmt) {
293 return 0;
294 // prb_ok=0;
295 } else if (fabs(psrch_id->xlngsnt - xlng[ysrt[0]]) > psrch_id->max_jmp) {
296 prb_ok=0;
297 } else {
298 psrch_id->call_code = -2;
299 return 1;
300 }
301 }
302 }
303
304 if (!prb_ok) {
305 if (psrch_id->brkset == 3) { //Try bracket method.
306 if (psrch_brckt(psrch_id)) {
307 psrch_id->call_code = -3;
308 return 1;
309
310 } else { // If bracket not ok, try downhill method.
311 psrch_dwnhll(psrch_id);
312 psrch_id->call_code = -4;
313 return 1;
314 }
315
316 } else { // Try downhill method if brackets not set.
317 psrch_dwnhll(psrch_id);
318 psrch_id->call_code = -4;
319 return 1;
320 }
321 }
322
323 return 1;
324 }
325
326 int psrch_getpt(psrch_vars *psrch_id, double *pt, double *yr)
327 /*
328 Returns 1 if search done, else 0.
329 */
330 {
331 int i=0,*ysrt,units;
332 double dflt,*x,*xlng,*y,strtpt;
333
334 // Speedup
335 dflt=psrch_id->dflt;
336 units=psrch_id->units;
337 strtpt=psrch_id->strtpt;
338 x=psrch_id->x;
339 xlng=psrch_id->xlng;
340 y=psrch_id->y;
341 ysrt=psrch_id->ysrt;
342
343 if (psrch_id->call_code != -1) { // Make sure the returned pt is the same is sent.
344 if (*pt != psrch_id->ptsnt) {
345 printf("%s%f%s%f%s%d\n",
346 "psrch_getpt error - point returned not same as sent: pt= ", *pt,
347 " != ptsnt= ", psrch_id->ptsnt,
348 " call_code= ", psrch_id->call_code);
349
350 // Set return variables from best available.
351 *pt=x[ysrt[0]];
352
353 // Make MAXIMA a MINIMA search.
354 if (psrch_id->srch_typ == MAXIMA) *yr=-psrch_id->y[ysrt[0]];
355 else *yr=psrch_id->y[ysrt[0]];
356
357 return 1;
358 }
359 }
360
361 if (psrch_id->srch_typ == MAXIMA) *yr = - *yr; // Make MAXIMA a MINIMA search.
362
363 if (psrch_id->call_code < -1) { // Usual minimization routine after initialization.
364
365 psrch_cbrckt(psrch_id, *yr, psrch_id->xlngsnt); // do bracket check
366
367 if (*yr < y[ysrt[2]]) { // If better than worst pt of three, add it. Use parabola.
368
369 printf("%s%f%s%f\n", "better pt added yr=", *yr, " pt= ", *pt);
370
371 y[ysrt[2]]=*yr;
372 xlng[ysrt[2]]=psrch_id->xlngsnt;
373 x[ysrt[2]]=*pt;
374
375 if (*yr < y[ysrt[0]]) { // Apply search gain measure
376 psrch_id->ygain = 1;
377 if ((y[ysrt[0]] - *yr) < psrch_id->srch_lmt)
378 psrch_id->imprv_fctr += psrch_id->prtl_imprv;
379 else psrch_id->imprv_fctr = 0.0;
380 } else psrch_id->imprv_fctr += psrch_id->prtl_imprv;
381
382 psrch_srt(psrch_id);
383
384 if (psrch_id->imprv_fctr > psrch_id->imprv_fctr_lmt) { // Search done.
385 *pt=x[ysrt[0]];
386 if (psrch_id->srch_typ == MAXIMA) *yr=-psrch_id->y[ysrt[0]];
387 else *yr=psrch_id->y[ysrt[0]];
388 printf("%s%f%s%f\n", "psrch_id->imprv_fctr > psrch_id->imprv_fctr_lmt yr=", *yr, " pt= ", *pt);
389 return 1;
390 }
391
392 if (!psrch_parest(psrch_id)) { // Search done.
393 *pt=x[ysrt[0]];
394 if (psrch_id->srch_typ == MAXIMA) *yr=-psrch_id->y[ysrt[0]];
395 else *yr=psrch_id->y[ysrt[0]];
396 return 1;
397 }
398
399 *pt=psrch_id->ptsnt=strtpt+psrch_id->xlngsnt;
400 return 0; // Get another yr from function searched.
401
402 } else if (psrch_id->brkset == 3) {
403 // Use Bracket or Downhill methods until better pt found
404 psrch_id->imprv_fctr
405 += (psrch_id->call_code == -2 ? 1.0 : psrch_id->prtl_imprv);
406 if (!psrch_brckt(psrch_id)) { // Search done
407 *pt=x[ysrt[0]];
408 if (psrch_id->srch_typ == MAXIMA) *yr=-psrch_id->y[ysrt[0]];
409 else *yr=psrch_id->y[ysrt[0]];
410 return 1;
411 }
412 *pt=psrch_id->ptsnt=strtpt+psrch_id->xlngsnt;
413
414 psrch_id->call_code = -3;
415 return 0; // Get another yr from function searched.
416
417 } else { // Downhill method
418 psrch_id->imprv_fctr
419 += (psrch_id->call_code == -2 ? 1.0 : psrch_id->prtl_imprv);
420 psrch_dwnhll(psrch_id);
421 *pt=psrch_id->ptsnt=strtpt+psrch_id->xlngsnt;
422
423 psrch_id->call_code = -4;
424 return 0; // Get another yr from function searched.
425 }
426
427 } else if (psrch_id->call_code == -1) { // Store and give pt back for