xue@11
|
1 /*
|
xue@11
|
2 Harmonic sinusoidal modelling and tools
|
xue@11
|
3
|
xue@11
|
4 C++ code package for harmonic sinusoidal modelling and relevant signal processing.
|
xue@11
|
5 Centre for Digital Music, Queen Mary, University of London.
|
xue@11
|
6 This file copyright 2011 Wen Xue.
|
xue@11
|
7
|
xue@11
|
8 This program is free software; you can redistribute it and/or
|
xue@11
|
9 modify it under the terms of the GNU General Public License as
|
xue@11
|
10 published by the Free Software Foundation; either version 2 of the
|
xue@11
|
11 License, or (at your option) any later version.
|
xue@11
|
12 */
|
xue@1
|
13 //---------------------------------------------------------------------------
|
xue@1
|
14
|
xue@1
|
15 #include <math.h>
|
xue@1
|
16 #include <memory.h>
|
xue@1
|
17 #include <stdlib.h>
|
xue@1
|
18 #include "opt.h"
|
xue@1
|
19 #include "matrix.h"
|
xue@1
|
20
|
Chris@5
|
21 /** @file opt.h - optimization routines */
|
Chris@5
|
22
|
xue@1
|
23 //---------------------------------------------------------------------------
|
xue@1
|
24 //macro nsign: judges if two double-precision floating pointer numbers have different signs
|
xue@1
|
25 #define nsign(x1, x2) (0x80&(((char*)&x1)[7]^((char*)&x2)[7]))
|
xue@1
|
26
|
xue@1
|
27 //---------------------------------------------------------------------------
|
Chris@5
|
28 /**
|
xue@1
|
29 function GradientMax: gradient method maximizing F(x)
|
xue@1
|
30
|
xue@1
|
31 In: function F(x) and its derivative dF(x)
|
xue@1
|
32 start[dim]: initial value of x
|
xue@1
|
33 ep: a stopping threshold
|
xue@1
|
34 maxiter: maximal number of iterates
|
xue@1
|
35 Out: start[dim]: the final x
|
xue@1
|
36
|
xue@1
|
37 Returnx max F(final x)
|
xue@1
|
38 */
|
xue@1
|
39 double GradientMax(void* params, double (*F)(int, double*, void*), void (*dF)(double*, int, double*, void*),
|
xue@1
|
40 int dim, double* start, double ep, int maxiter)
|
xue@1
|
41 {
|
xue@1
|
42 double oldG, G=F(dim, start, params), *dG=new double[dim], *X=new double[dim];
|
xue@1
|
43
|
xue@1
|
44 int iter=0;
|
xue@1
|
45 do
|
xue@1
|
46 {
|
xue@1
|
47 oldG=G;
|
xue@1
|
48
|
xue@1
|
49 dF(dG, dim, start, params);
|
xue@1
|
50 Search1Dmaxfrom(params, F, dim, start, dG, 1e-2, 1e-6, 1e6, X, &G, ep);
|
xue@1
|
51 memcpy(start, X, sizeof(double)*dim);
|
xue@1
|
52
|
xue@1
|
53 iter++;
|
xue@1
|
54 }
|
xue@1
|
55 while (iter<maxiter && G>oldG*(1+ep));
|
xue@1
|
56 return G;
|
xue@1
|
57 }//GradientMax
|
xue@1
|
58
|
Chris@5
|
59 /**
|
xue@1
|
60 function GradientOneMax: gradient method maximizing F(x), which searches $dim dimensions one after
|
xue@1
|
61 another instead of taking the gradient descent direction.
|
xue@1
|
62
|
xue@1
|
63 In: function F(x) and its derivative dF(x)
|
xue@1
|
64 start[dim]: initial value of x
|
xue@1
|
65 ep: a stopping threshold
|
xue@1
|
66 maxiter: maximal number of iterates
|
xue@1
|
67 Out: start[dim]: the final x
|
xue@1
|
68
|
xue@1
|
69 Returnx max F(final x)
|
xue@1
|
70 */
|
xue@1
|
71 double GradientOneMax(void* params, double (*F)(int, double*, void*), void (*dF)(double*, int, double*, void*),
|
xue@1
|
72 int dim, double* start, double ep, int maxiter)
|
xue@1
|
73 {
|
xue@1
|
74 double oldG, G=F(dim, start, params), *dG=new double[dim], *X=new double[dim];
|
xue@1
|
75
|
xue@1
|
76 int iter=0;
|
xue@1
|
77 do
|
xue@1
|
78 {
|
xue@1
|
79 oldG=G;
|
xue@1
|
80 for (int d=0; d<dim; d++)
|
xue@1
|
81 {
|
xue@1
|
82 dF(dG, dim, start, params);
|
xue@1
|
83 for (int i=0; i<dim; i++) if (d!=i) dG[i]=0;
|
xue@1
|
84 if (dG[d]!=0)
|
xue@1
|
85 {
|
xue@1
|
86 Search1Dmaxfrom(params, F, dim, start, dG, 1e-2, 1e-6, 1e6, X, &G, ep);
|
xue@1
|
87 memcpy(start, X, sizeof(double)*dim);
|
xue@1
|
88 }
|
xue@1
|
89 }
|
xue@1
|
90 iter++;
|
xue@1
|
91 }
|
xue@1
|
92 while (iter<maxiter && G>oldG*(1+ep));
|
xue@1
|
93 return G;
|
xue@1
|
94 }//GradientOneMax
|
xue@1
|
95
|
xue@1
|
96 //---------------------------------------------------------------------------
|
Chris@5
|
97 /**
|
xue@1
|
98 function mindenormaldouble: returns the minimal denormal double-precision floating point number.
|
xue@1
|
99 */
|
xue@1
|
100 double mindenormaldouble()
|
xue@1
|
101 {
|
xue@1
|
102 double result=0;
|
xue@1
|
103 *((unsigned char*)&result)=1;
|
xue@1
|
104 return result;
|
xue@1
|
105 }//mindenormaldouble
|
xue@1
|
106
|
Chris@5
|
107 /**
|
xue@1
|
108 function minnormaldouble: returns the minimal normal double-precision floating point number.
|
xue@1
|
109 */
|
xue@1
|
110 double minnormaldouble()
|
xue@1
|
111 {
|
xue@1
|
112 double result=0;
|
xue@1
|
113 ((unsigned char*)&result)[6]=16;
|
xue@1
|
114 return result;
|
xue@1
|
115 }//minnormaldouble
|
xue@1
|
116
|
xue@1
|
117 /*
|
xue@1
|
118 //function nextdouble: returns the next double value after x in the direction of dir
|
xue@1
|
119 */
|
xue@1
|
120 double nextdouble(double x, double dir)
|
xue@1
|
121 {
|
xue@1
|
122 if (x==0)
|
xue@1
|
123 {
|
xue@1
|
124 if (dir<0) return -mindenormaldouble();
|
xue@1
|
125 else return mindenormaldouble();
|
xue@1
|
126 }
|
xue@1
|
127 else if ((x>0)^(dir>0)) *((__int64*)&x)-=1;
|
xue@1
|
128 else *((__int64*)&x)+=1;
|
xue@1
|
129 return x;
|
xue@1
|
130 }//nextdouble
|
xue@1
|
131
|
Chris@5
|
132 /**
|
xue@1
|
133 function Newton: Newton method for finding the root of y(x)
|
xue@1
|
134
|
xue@1
|
135 In: function y(x) and its derivative y'(x)
|
xue@1
|
136 params: additional argument for calling y and y'
|
xue@1
|
137 x0: inital value of x
|
xue@1
|
138 epx, epy: stopping thresholds on x and y respectively
|
xue@1
|
139 maxiter: maximal number of iterates
|
xue@1
|
140 Out: x0: the root
|
xue@1
|
141
|
xue@1
|
142 Returns 0 if successful.
|
xue@1
|
143 */
|
xue@1
|
144 double Newton(double& x0, double (*y)(double, void*), double (*y1)(double, void*), void* params, double epx, int maxiter, double epy)
|
xue@1
|
145 {
|
xue@1
|
146 double y0=y(x0, params), x2, y2, dx, dy, x3, y3;
|
xue@1
|
147 if (y0==0) return 0;
|
xue@1
|
148
|
xue@1
|
149 int iter=0;
|
xue@1
|
150 while (iter<maxiter)
|
xue@1
|
151 {
|
xue@1
|
152 dy=y1(x0, params);
|
xue@1
|
153 if (dy==0) return -1;
|
xue@1
|
154 x2=x0-y0/dy;
|
xue@1
|
155 y2=y(x2, params);
|
xue@1
|
156 while (fabs(y2)>fabs(y0))
|
xue@1
|
157 {
|
xue@1
|
158 x2=(x2+x0)/2;
|
xue@1
|
159 y2=y(x2, params);
|
xue@1
|
160 }
|
xue@1
|
161 if (y2==0)
|
xue@1
|
162 {
|
xue@1
|
163 x0=x2;
|
xue@1
|
164 return 0;
|
xue@1
|
165 }
|
xue@1
|
166 dx=fabs(x0-x2);
|
xue@1
|
167 if (dx==0)
|
xue@1
|
168 {
|
xue@1
|
169 if (fabs(y0)<epy) return 0;
|
xue@1
|
170 else return -1;
|
xue@1
|
171 }
|
xue@1
|
172 else if (dx<epx)
|
xue@1
|
173 {
|
xue@1
|
174 if nsign(y2, y0)
|
xue@1
|
175 {
|
xue@1
|
176 x0=x2;
|
xue@1
|
177 return dx;
|
xue@1
|
178 }
|
xue@1
|
179 else
|
xue@1
|
180 {
|
xue@1
|
181 x3=2*x2-x0;
|
xue@1
|
182 y3=y(x3, params);
|
xue@1
|
183 if (y3==0)
|
xue@1
|
184 {
|
xue@1
|
185 x0=x3;
|
xue@1
|
186 return 0;
|
xue@1
|
187 }
|
xue@1
|
188 if nsign(y2, y3)
|
xue@1
|
189 {
|
xue@1
|
190 x0=x2;
|
xue@1
|
191 return dx;
|
xue@1
|
192 }
|
xue@1
|
193 else if (fabs(y3)<fabs(y2))
|
xue@1
|
194 {
|
xue@1
|
195 x2=x3, y2=y3;
|
xue@1
|
196 }
|
xue@1
|
197 }
|
xue@1
|
198 }
|
xue@1
|
199 x0=x2;
|
xue@1
|
200 y0=y2;
|
xue@1
|
201 iter++;
|
xue@1
|
202 }
|
xue@1
|
203 if (fabs(y0)<epy) return 0;
|
xue@1
|
204 else return -1;
|
xue@1
|
205 }//Newton
|
xue@1
|
206
|
Chris@5
|
207 /**
|
xue@1
|
208 function Newton: Newton method for finding the root of y(x), for use when it's more efficient for y'
|
xue@1
|
209 and y to be calculated in one function call
|
xue@1
|
210
|
xue@1
|
211 In: function dy(x) that computes y and y' at x
|
xue@1
|
212 params: additional argument for calling dy
|
xue@1
|
213 yshift: the byte shift in params where dy(x) put the value of y(x) as double type
|
xue@1
|
214 x0: inital value of x
|
xue@1
|
215 epx, epy: stopping thresholds on x and y respectively
|
xue@1
|
216 maxiter: maximal number of iterates
|
xue@1
|
217 Out: x0: the root
|
xue@1
|
218
|
xue@1
|
219 Returns 0 if successful.
|
xue@1
|
220 */
|
xue@1
|
221 int Newton(double& x0, double (*dy)(double, void*), void* params, int yshift, double epx, int maxiter, double epy, double xmin, double xmax)
|
xue@1
|
222 {
|
xue@1
|
223 double y0, x2, y2, dx, dy0, dy2, dy3, x3, y3;
|
xue@1
|
224 dy0=dy(x0, params);
|
xue@1
|
225 y0=*(double*)((char*)params+yshift);
|
xue@1
|
226 if (y0==0) return 0;
|
xue@1
|
227
|
xue@1
|
228 int iter=0;
|
xue@1
|
229 while (iter<maxiter)
|
xue@1
|
230 {
|
xue@1
|
231 if (dy0==0) return -1;
|
xue@1
|
232 dx=-y0/dy0;
|
xue@1
|
233
|
xue@1
|
234 x2=x0+dx;
|
xue@1
|
235 dy2=dy(x2, params);
|
xue@1
|
236 y2=*(double*)((char*)params+yshift);
|
xue@1
|
237 //make sure the iteration yields a y closer to 0
|
xue@1
|
238 while (fabs(y2)>fabs(y0))
|
xue@1
|
239 {
|
xue@1
|
240 dx/=2;
|
xue@1
|
241 x2=x0+dx;
|
xue@1
|
242 dy2=dy(x2, params);
|
xue@1
|
243 y2=*(double*)((char*)params+yshift);
|
xue@1
|
244 }
|
xue@1
|
245 //the above may fail due to precision problem, i.e. x0+dx=x0
|
xue@1
|
246 if (x2==x0)
|
xue@1
|
247 {
|
xue@1
|
248 if (fabs(y0)<epy) return 0;
|
xue@1
|
249 else
|
xue@1
|
250 {
|
xue@1
|
251 x2=nextdouble(x0, dx);
|
xue@1
|
252 dy2=dy(x2, params);
|
xue@1
|
253 y2=*(double*)((char*)params+yshift);
|
xue@1
|
254 if (y2==0)
|
xue@1
|
255 {
|
xue@1
|
256 x0=x2;
|
xue@1
|
257 return 0;
|
xue@1
|
258 }
|
xue@1
|
259 else
|
xue@1
|
260 {
|
xue@1
|
261 dx=x2-x0;
|
xue@1
|
262 while (y2==y0)
|
xue@1
|
263 {
|
xue@1
|
264 dx*=2;
|
xue@1
|
265 x2=x0+dx;
|
xue@1
|
266 dy2=dy(x2, params);
|
xue@1
|
267 y2=*(double*)((char*)params+yshift);
|
xue@1
|
268 }
|
xue@1
|
269 if nsign(y2, y0)
|
xue@1
|
270 return fabs(dx);
|
xue@1
|
271 else if (fabs(y2)>fabs(y0))
|
xue@1
|
272 return -1;
|
xue@1
|
273 }
|
xue@1
|
274 }
|
xue@1
|
275 }
|
xue@1
|
276 else if (y2==0)
|
xue@1
|
277 {
|
xue@1
|
278 x0=x2;
|
xue@1
|
279 return 0;
|
xue@1
|
280 }
|
xue@1
|
281 else
|
xue@1
|
282 {
|
xue@1
|
283 double fabsdx=fabs(dx);
|
xue@1
|
284 if (fabsdx<epx)
|
xue@1
|
285 {
|
xue@1
|
286 if nsign(y2, y0)
|
xue@1
|
287 {
|
xue@1
|
288 x0=x2;
|
xue@1
|
289 return fabsdx;
|
xue@1
|
290 }
|
xue@1
|
291 else
|
xue@1
|
292 {
|
xue@1
|
293 x3=x2+dx;
|
xue@1
|
294 dy3=dy(x3, params);
|
xue@1
|
295 y3=*(double*)((char*)params+yshift);
|
xue@1
|
296 if (y3==0){x0=x3; return 0;}
|
xue@1
|
297 if nsign(y2, y3)
|
xue@1
|
298 {
|
xue@1
|
299 x0=x2;
|
xue@1
|
300 return dx;
|
xue@1
|
301 }
|
xue@1
|
302 else if (fabs(y3)<fabs(y2))
|
xue@1
|
303 {
|
xue@1
|
304 x2=x3, y2=y3, dy2=dy3;
|
xue@1
|
305 }
|
xue@1
|
306 }
|
xue@1
|
307 }
|
xue@1
|
308 }
|
xue@1
|
309 if (x2<xmin || x2>xmax) return -1;
|
xue@1
|
310 x0=x2;
|
xue@1
|
311 y0=y2;
|
xue@1
|
312 dy0=dy2;
|
xue@1
|
313 iter++;
|
xue@1
|
314 }
|
xue@1
|
315 if (fabs(y0)<epy) return 0;
|
xue@1
|
316 else return -1;
|
xue@1
|
317 }//Newton
|
xue@1
|
318
|
Chris@5
|
319 /**
|
xue@1
|
320 function init_Newton_1d: finds an initial x and interval (xa, xb) from which to search for an extremum
|
xue@1
|
321
|
xue@1
|
322 On calling, xa<x0<xb, y(x0)>=y(xa), y(x0)>=y(xb) (or y(x0)<=y(xa), y(x0)<=y(xb)).
|
xue@1
|
323 On return, y(x0)>=y(xa), y(x0)>=y(xb) (or y(x0)<=y(xa), y(x0)<=y(xb)), and either y'(xa)>0, y'(xb)<0
|
xue@1
|
324 (or y'(xa)<0, y'(xb)>0), or xb-xa<epx.
|
xue@1
|
325
|
xue@1
|
326 In: function dy(x) that computes y and y' at x
|
xue@1
|
327 params: additional argument for calling dy
|
xue@1
|
328 yshift: the byte shift in params where dy(x) put the value of y(x) as double type
|
xue@1
|
329 x0: inital value of x
|
xue@1
|
330 (xa, xb): initial search range
|
xue@1
|
331 epx: tolerable error of extremum
|
xue@1
|
332 Out: improved initial value x0 and search range (xa, xb)
|
xue@1
|
333
|
xue@1
|
334 Returns xb-xa if successful, -0.5 or -0.6 if y(x0) is not bigger(smaller) than both y(xa) and y(xb)
|
xue@1
|
335 */
|
xue@1
|
336 double init_Newton_1d(double& x0, double& xa, double& xb, double (*dy)(double, void*), void* params, int yshift, double epx)
|
xue@1
|
337 {
|
xue@1
|
338 double dy0, dya, dyb, y0, ya, yb;
|
xue@1
|
339 dy0=dy(x0, params); y0=*(double*)((char*)params+yshift);
|
xue@1
|
340 dya=dy(xa, params); ya=*(double*)((char*)params+yshift);
|
xue@1
|
341 dyb=dy(xb, params); yb=*(double*)((char*)params+yshift);
|
xue@1
|
342 if (y0>=ya && y0>=yb)
|
xue@1
|
343 {
|
xue@1
|
344 //look for xa<x0<xb, ya<y0>yb, so that dya>0, dyb<0
|
xue@1
|
345 while ((dya<=0 || dyb>=0) && xb-xa>=epx)
|
xue@1
|
346 {
|
xue@1
|
347 double x1=(x0+xa)/2;
|
xue@1
|
348 double dy1=dy(x1, params);
|
xue@1
|
349 double y1=*(double*)((char*)params+yshift);
|
xue@1
|
350 double x2=(x0+xb)/2;
|
xue@1
|
351 double dy2=dy(x2, params);
|
xue@1
|
352 double y2=*(double*)((char*)params+yshift);
|
xue@1
|
353 if (y0>=y1 && y0>=y2) xa=x1, ya=y1, dya=dy1, xb=x2, yb=y2, dyb=dy2;
|
xue@1
|
354 else if (y1>=y0 && y1>=y2) {xb=x0, yb=y0, dyb=dy0; x0=x1, y0=y1, dy0=dy1;}
|
xue@1
|
355 else {xa=x0, ya=y0, dya=dy0; x0=x2, y0=y2, dy0=dy2;}
|
xue@1
|
356 }
|
xue@1
|
357 return xb-xa;
|
xue@1
|
358 }
|
xue@1
|
359 else if (y0<=ya && y0<=yb)
|
xue@1
|
360 {
|
xue@1
|
361 //look for xa<x0<xb, ya>y0<yb, so that dya<0, dyb>0
|
xue@1
|
362 while ((dya>=0 || dyb<=0) && xb-xa>=epx)
|
xue@1
|
363 {
|
xue@1
|
364 double x1=(x0+xa)/2;
|
xue@1
|
365 double dy1=dy(x1, params);
|
xue@1
|
366 double y1=*(double*)((char*)params+yshift);
|
xue@1
|
367 double x2=(x0+xb)/2;
|
xue@1
|
368 double dy2=dy(x2, params);
|
xue@1
|
369 double y2=*(double*)((char*)params+yshift);
|
xue@1
|
370 if (y0<=y1 && y0<=y2) xa=x1, ya=y1, dya=dy1, xb=x2, yb=y2, dyb=dy2;
|
xue@1
|
371 else if (y1<=y0 && y1<=y2) {xb=x0, yb=y0, dyb=dy0; x0=x1, y0=y1, dy0=dy1;}
|
xue@1
|
372 else {xa=x0, ya=y0, dya=dy0; x0=x2, y0=y2, dy0=dy2;}
|
xue@1
|
373 }
|
xue@1
|
374 return xb-xa;
|
xue@1
|
375 }
|
xue@1
|
376 else//meaning y(x0) is not bigger(smaller) than both y(xa) and y(xb)
|
xue@1
|
377 {
|
xue@1
|
378 if (y0<=yb) return -0.5; //ya<=y0<=yb
|
xue@1
|
379 if (y0<=ya) return -0.6; //ya>=y0>=yb
|
xue@1
|
380 }
|
xue@1
|
381 return -0;
|
xue@1
|
382 }//init_Newton_1d
|
xue@1
|
383
|
Chris@5
|
384 /**
|
xue@1
|
385 function init_Newton_1d_max: finds an initial x and interval (xa, xb) from which to search for an
|
xue@1
|
386 maximum.
|
xue@1
|
387
|
xue@1
|
388 On calling xa<x0<xb, and it is assumed that y(x0)>=y(xa), y(x0)>=y(xb). If the conditions on y are met,
|
xue@1
|
389 it find a new set of xa<x0<xb, so that y(x0)>y(xa), y(x0)>=y(xb), y'(xa)>0, y'(xb)<0, or xb-xa<epx,
|
xue@1
|
390 and returns xb-xa>0; otherwise it returns a negative value showing the initial order of y(x0), y(xa)
|
xue@1
|
391 and y(xb) as
|
xue@1
|
392 -0.5 for ya<=y0<yb
|
xue@1
|
393 -0.6 for yb<=y0<ya
|
xue@1
|
394 -0.7 for y0<ya<yb
|
xue@1
|
395 -0.8 for y0<yb<=ya
|
xue@1
|
396
|
xue@1
|
397 In: function dy(x) that computes y and y' at x
|
xue@1
|
398 params: additional argument for calling dy
|
xue@1
|
399 yshift: the byte shift in params where dy(x) put the value of y(x) as double type
|
xue@1
|
400 x0: inital value of x
|
xue@1
|
401 (xa, xb): initial search range
|
xue@1
|
402 epx: tolerable error of extremum
|
xue@1
|
403 Out: improved initial value x0 and search range (xa, xb)
|
xue@1
|
404
|
xue@1
|
405 Returns xb-xa if successful, a negative value as above if not.
|
xue@1
|
406 */
|
xue@1
|
407 double init_Newton_1d_max(double& x0, double& xa, double& xb, double (*dy)(double, void*), void* params, int yshift, double epx)
|
xue@1
|
408 {
|
xue@1
|
409 double dy0, dya, dyb, y0, ya, yb;
|
xue@1
|
410 dy0=dy(x0, params); y0=*(double*)((char*)params+yshift);
|
xue@1
|
411 dya=dy(xa, params); ya=*(double*)((char*)params+yshift);
|
xue@1
|
412 dyb=dy(xb, params); yb=*(double*)((char*)params+yshift);
|
xue@1
|
413 if (y0>=ya && y0>=yb)
|
xue@1
|
414 {
|
xue@1
|
415 //look for xa<x0<xb, ya<y0>yb, so that dya>0, dyb<0
|
xue@1
|
416 while ((dya<=0 || dyb>=0) && xb-xa>=epx)
|
xue@1
|
417 {
|
xue@1
|
418 double x1=(x0+xa)/2;
|
xue@1
|
419 double dy1=dy(x1, params);
|
xue@1
|
420 double y1=*(double*)((char*)params+yshift);
|
xue@1
|
421 double x2=(x0+xb)/2;
|
xue@1
|
422 double dy2=dy(x2, params);
|
xue@1
|
423 double y2=*(double*)((char*)params+yshift);
|
xue@1
|
424 if (y0>=y1 && y0>=y2) xa=x1, ya=y1, dya=dy1, xb=x2, yb=y2, dyb=dy2;
|
xue@1
|
425 else if (y1>=y0 && y1>=y2) {xb=x0, yb=y0, dyb=dy0; x0=x1, y0=y1, dy0=dy1;}
|
xue@1
|
426 else {xa=x0, ya=y0, dya=dy0; x0=x2, y0=y2, dy0=dy2;}
|
xue@1
|
427 }
|
xue@1
|
428 return xb-xa;
|
xue@1
|
429 }
|
xue@1
|
430 else//meaning y(x0) is not bigger than both y(xa) and y(xb)
|
xue@1
|
431 {
|
xue@1
|
432 if (y0<yb && y0>=ya) return -0.5; //ya<=y0<yb
|
xue@1
|
433 if (y0<ya && y0>=yb) return -0.6; //ya>y0>=yb
|
xue@1
|
434 if (ya<yb) return -0.7; //y0<ya<yb
|
xue@1
|
435 return -0.8; //y0<yb<=ya
|
xue@1
|
436 }
|
xue@1
|
437 }//init_Newton_1d_max
|
xue@1
|
438
|
Chris@5
|
439 /**
|
xue@1
|
440 function init_Newton_1d_min: finds an initial x and interval (xa, xb) from which to search for an
|
xue@1
|
441 manimum.
|
xue@1
|
442
|
xue@1
|
443 On calling xa<x0<xb, and it is assumed that y(x0)<=y(xa), y(x0)<=y(xb). If the conditions on y are met,
|
xue@1
|
444 it find a new set of xa<x0<xb, so that y(x0)<y(xa), y(x0)<=y(xb), y'(xa)<0, y'(xb)>0, or xb-xa<epx,
|
xue@1
|
445 and returns xb-xa>0; otherwise it returns a negative value showing the initial order of y(x0), y(xa)
|
xue@1
|
446 and y(xb) as
|
xue@1
|
447 -0.5 for ya<=y0<yb
|
xue@1
|
448 -0.6 for yb<=y0<ya
|
xue@1
|
449 -0.7 for y0>yb>ya
|
xue@1
|
450 -0.8 for y0>ya>=yb
|
xue@1
|
451
|
xue@1
|
452 In: function dy(x) that computes y and y' at x
|
xue@1
|
453 params: additional argument for calling dy
|
xue@1
|
454 yshift: the byte shift in params where dy(x) put the value of y(x) as double type
|
xue@1
|
455 x0: inital value of x
|
xue@1
|
456 (xa, xb): initial search range
|
xue@1
|
457 epx: tolerable error of extremum
|
xue@1
|
458 Out: improved initial value x0 and search range (xa, xb)
|
xue@1
|
459
|
xue@1
|
460 Returns xb-xa if successful, a negative value as above if not.
|
xue@1
|
461 */
|
xue@1
|
462 double init_Newton_1d_min(double& x0, double& xa, double& xb, double (*dy)(double, void*), void* params, int yshift, double epx)
|
xue@1
|
463 {
|
xue@1
|
464 double dy0, dya, dyb, y0, ya, yb;
|
xue@1
|
465 dy0=dy(x0, params); y0=*(double*)((char*)params+yshift);
|
xue@1
|
466 dya=dy(xa, params); ya=*(double*)((char*)params+yshift);
|
xue@1
|
467 dyb=dy(xb, params); yb=*(double*)((char*)params+yshift);
|
xue@1
|
468 if (y0<=ya && y0<=yb)
|
xue@1
|
469 {
|
xue@1
|
470 //look for xa<x0<xb, ya>y0<yb, so that dya<0, dyb>0
|
xue@1
|
471 while ((dya>=0 || dyb<=0) && xb-xa>=epx)
|
xue@1
|
472 {
|
xue@1
|
473 double x1=(x0+xa)/2;
|
xue@1
|
474 double dy1=dy(x1, params);
|
xue@1
|
475 double y1=*(double*)((char*)params+yshift);
|
xue@1
|
476 double x2=(x0+xb)/2;
|
xue@1
|
477 double dy2=dy(x2, params);
|
xue@1
|
478 double y2=*(double*)((char*)params+yshift);
|
xue@1
|
479 if (y0<=y1 && y0<=y2) xa=x1, ya=y1, dya=dy1, xb=x2, yb=y2, dyb=dy2;
|
xue@1
|
480 else if (y1<=y0 && y1<=y2) {xb=x0, yb=y0, dyb=dy0; x0=x1, y0=y1, dy0=dy1;}
|
xue@1
|
481 else {xa=x0, ya=y0, dya=dy0; x0=x2, y0=y2, dy0=dy2;}
|
xue@1
|
482 }
|
xue@1
|
483 return xb-xa;
|
xue@1
|
484 }
|
xue@1
|
485 else//meaning y(x0) is not bigger than both y(xa) and y(xb)
|
xue@1
|
486 {
|
xue@1
|
487 if (y0<yb && y0>=ya) return -0.5; //ya<=y0<yb
|
xue@1
|
488 if (y0<ya && y0>=yb) return -0.6; //ya>y0>=yb
|
xue@1
|
489 if (ya<yb) return -0.7; //ya<yb<y0
|
xue@1
|
490 return -0.8; //yb<=ya<y0
|
xue@1
|
491 }
|
xue@1
|
492 }//init_Newton_1d_min
|
xue@1
|
493
|
Chris@5
|
494 /**
|
xue@1
|
495 function Newton_1d_max: Newton method for maximizing y(x). On calling y(x0)>=y(xa), y(x0)>=y(xb),
|
xue@1
|
496 y'(xa)>0, y'(xb)<0
|
xue@1
|
497
|
xue@1
|
498 In: function ddy(x) that computes y, y' and y" at x
|
xue@1
|
499 params: additional argument for calling ddy
|
xue@1
|
500 y1shift: the byte shift in params where ddy(x) put the value of y'(x) as double type
|
xue@1
|
501 yshift: the byte shift in params where ddy(x) put the value of y(x) as double type
|
xue@1
|
502 x0: inital value of x
|
xue@1
|
503 (xa, xb): initial search range
|
xue@1
|
504 epx, epdy: stopping threshold for x and y', respectively
|
xue@1
|
505 maxiter: maximal number of iterates
|
xue@1
|
506 Out: x0: the maximum
|
xue@1
|
507
|
xue@1
|
508 Returns the error bound for the maximum x0, or
|
xue@1
|
509 -2 if y(*) do not satisfy input conditions
|
xue@1
|
510 -3 if y'(*) do not satisfy input conditions
|
xue@1
|
511 */
|
xue@1
|
512 double Newton_1d_max(double& x0, double xa, double xb, double (*ddy)(double, void*), void* params, int y1shift, int yshift, double epx, int maxiter, double epdy, bool checkinput)
|
xue@1
|
513 {
|
xue@1
|
514 double x1, y1, x2, y2, dx, dy1, dy2, dy3, x3, y3, ddy1, ddy2, ddy3;
|
xue@1
|
515 double y0, ya, yb, dy0, dya, dyb, ddy0;
|
xue@1
|
516 if (xb-xa<epx) {return xb-xa;}
|
xue@1
|
517 ddy(xa, params);
|
xue@1
|
518 dya=*(double*)((char*)params+y1shift), ya=*(double*)((char*)params+yshift);
|
xue@1
|
519 ddy(xb, params);
|
xue@1
|
520 dyb=*(double*)((char*)params+y1shift), yb=*(double*)((char*)params+yshift);
|
xue@1
|
521 ddy0=ddy(x0, params);
|
xue@1
|
522 dy0=*(double*)((char*)params+y1shift), y0=*(double*)((char*)params+yshift);
|
xue@1
|
523
|
xue@1
|
524 if (checkinput)
|
xue@1
|
525 {
|
xue@1
|
526 if (y0>=ya && y0>=yb)
|
xue@1
|
527 {
|
xue@1
|
528 if (dya<=0 || dyb>=0) return -3;
|
xue@1
|
529 }
|
xue@1
|
530 else return -2;
|
xue@1
|
531 }
|
xue@1
|
532
|
xue@1
|
533 if (dy0==0) return 0;
|
xue@1
|
534
|
xue@1
|
535 int iter=0;
|
xue@1
|
536 while (iter<maxiter)
|
xue@1
|
537 {
|
xue@1
|
538 if (ddy0>=0) //newton method not applicable
|
xue@1
|
539 {
|
xue@1
|
540 do
|
xue@1
|
541 {
|
xue@1
|
542 x1=(x0+xa)/2;
|
xue@1
|
543 ddy1=ddy(x1, params);
|
xue@1
|
544 dy1=*(double*)((char*)params+y1shift), y1=*(double*)((char*)params+yshift);
|
xue@1
|
545 x2=(x0+xb)/2;
|
xue@1
|
546 ddy2=ddy(x2, params);
|
xue@1
|
547 dy2=*(double*)((char*)params+y1shift), y2=*(double*)((char*)params+yshift);
|
xue@1
|
548 if (y0>=y1 && y0>=y2) xa=x1, ya=y1, dya=dy1, xb=x2, yb=y2, dyb=dy2;
|
xue@1
|
549 else if (y1>=y0 && y1>=y2) {xb=x0, yb=y0, dyb=dy0; x0=x1, y0=y1, dy0=dy1, ddy0=ddy1;}
|
xue@1
|
550 else {xa=x0, ya=y0, dya=dy0; x0=x2, y0=y2, dy0=dy2, ddy0=ddy2;}
|
xue@1
|
551 }
|
xue@1
|
552 while(ddy0>=0 && (dya<=0 || dyb>=0) && xb-xa>=epx);
|
xue@1
|
553 if (xb-xa<epx) return xb-xa;
|
xue@1
|
554 }
|
xue@1
|
555 if (fabs(dy0)<epdy) return 0;
|
xue@1
|
556 //x2: Newton point
|
xue@1
|
557 //dx is the same direction as dy0 if ddy0<0
|
xue@1
|
558 dx=-dy0/ddy0;
|
xue@1
|
559 x2=x0+dx;
|
xue@1
|
560 if (x2==x0) return 0;
|
xue@1
|
561
|
xue@1
|
562 //make sure x2 lies in [xa, xb]
|
xue@1
|
563 if (x2<xa) {x2=(xa+x0)/2; dx=x2-x0;}
|
xue@1
|
564 else if (x2>xb) {x2=(xb+x0)/2; dx=x2-x0;}
|
xue@1
|
565 ddy2=ddy(x2, params);
|
xue@1
|
566 y2=*(double*)((char*)params+yshift);
|
xue@1
|
567
|
xue@1
|
568 //make sure the iteration yields a larger y
|
xue@1
|
569 if (y2<y0)
|
xue@1
|
570 {
|
xue@1
|
571 while (y2<y0 && dx>epx)
|
xue@1
|
572 {
|
xue@1
|
573 dx/=2;
|
xue@1
|
574 x2=x0+dx;
|
xue@1
|
575 ddy2=ddy(x2, params);
|
xue@1
|
576 y2=*(double*)((char*)params+yshift);
|
xue@1
|
577 }
|
xue@1
|
578 if (y2<y0) return fabs(dx);
|
xue@1
|
579 }
|
xue@1
|
580
|
xue@1
|
581 //new x2 lies between xa and xb and y2>=y0
|
xue@1
|
582 dy2=*(double*)((char*)params+y1shift);
|
xue@1
|
583 if (fabs(dy2)<epdy){x0=x2; return 0;}
|
xue@1
|
584
|
xue@1
|
585 double fabsdx=fabs(dx);
|
xue@1
|
586 if (fabsdx<epx)
|
xue@1
|
587 {
|
xue@1
|
588 if nsign(dy0, dy2){x0=x2; return fabsdx;}
|
xue@1
|
589 x3=x2+dx;
|
xue@1
|
590 ddy3=ddy(x3, params);
|
xue@1
|
591 dy3=*(double*)((char*)params+y1shift), y3=*(double*)((char*)params+yshift);
|
xue@1
|
592 if (dy3==0) {x0=x3; return 0;}
|
xue@1
|
593 if (y2>=y3) {x0=x2; return fabsdx;}
|
xue@1
|
594 else x2=x3, y2=y3, dy2=dy3, ddy2=ddy3;
|
xue@1
|
595 }
|
xue@1
|
596
|
xue@1
|
597 if (x0<x2) xa=x0, ya=y0, dya=dy0;
|
xue@1
|
598 else xb=x0, yb=y0, dyb=dy0;
|
xue@1
|
599 x0=x2, y0=y2, dy0=dy2, ddy0=ddy2;
|
xue@1
|
600 iter++;
|
xue@1
|
601 }
|
xue@1
|
602 if (fabs(dy0)<epdy) return 0;
|
xue@1
|
603 else return -1;
|
xue@1
|
604 }//Newton_1d_max
|
xue@1
|
605
|
Chris@5
|
606 /**
|
xue@10
|
607 function Newton_1d_min: Newton method for minimizing y(x). On calling y(x0)<=y(xa), y(x0)<=y(xb),
|
xue@1
|
608 y'(xa)<0, y'(xb)>0
|
xue@1
|
609
|
xue@1
|
610 In: function ddy(x) that computes y, y' and y" at x
|
xue@1
|
611 params: additional argument for calling ddy
|
xue@1
|
612 y1shift: the byte shift in params where ddy(x) put the value of y'(x) as double type
|
xue@1
|
613 yshift: the byte shift in params where ddy(x) put the value of y(x) as double type
|
xue@1
|
614 x0: inital value of x
|
xue@1
|
615 (xa, xb): initial search range
|
xue@1
|
616 epx, epdy: stopping threshold for x and y', respectively
|
xue@1
|
617 maxiter: maximal number of iterates
|
xue@1
|
618 Out: x0: the maximum
|
xue@1
|
619
|
xue@1
|
620 Returns the error bound for the maximum x0, or
|
xue@1
|
621 -2 if y(*) do not satisfy input conditions
|
xue@1
|
622 -3 if y'(*) do not satisfy input conditions
|
xue@1
|
623 */
|
xue@1
|
624 double Newton_1d_min(double& x0, double xa, double xb, double (*ddy)(double, void*), void* params, int y1shift, int yshift, double epx, int maxiter, double epdy, bool checkinput)
|
xue@1
|
625 {
|
xue@1
|
626 double x1, y1, x2, y2, dx, dy1, dy2, dy3, x3, y3, ddy1, ddy2, ddy3;
|
xue@1
|
627 double y0, ya, yb, dy0, dya, dyb, ddy0;
|
xue@1
|
628 if (xb-xa<epx) {return xb-xa;}
|
xue@1
|
629 ddy(xa, params);
|
xue@1
|
630 dya=*(double*)((char*)params+y1shift), ya=*(double*)((char*)params+yshift);
|
xue@1
|
631 ddy(xb, params);
|
xue@1
|
632 dyb=*(double*)((char*)params+y1shift), yb=*(double*)((char*)params+yshift);
|
xue@1
|
633 ddy0=ddy(x0, params);
|
xue@1
|
634 dy0=*(double*)((char*)params+y1shift), y0=*(double*)((char*)params+yshift);
|
xue@1
|
635
|
xue@1
|
636 if (checkinput)
|
xue@1
|
637 {
|
xue@1
|
638 if (y0<=ya && y0<=yb)
|
xue@1
|
639 {
|
xue@1
|
640 if (dya>=0 || dyb<=0) return -3;
|
xue@1
|
641 }
|
xue@1
|
642 else return -2;
|
xue@1
|
643 }
|
xue@1
|
644
|
xue@1
|
645 if (dy0==0) return 0;
|
xue@1
|
646
|
xue@1
|
647 int iter=0;
|
xue@1
|
648 while (iter<maxiter)
|
xue@1
|
649 {
|
xue@1
|
650 if (ddy0<=0) //newton method not applicable
|
xue@1
|
651 {
|
xue@1
|
652 do
|
xue@1
|
653 {
|
xue@1
|
654 x1=(x0+xa)/2;
|
xue@1
|
655 ddy1=ddy(x1, params);
|
xue@1
|
656 dy1=*(double*)((char*)params+y1shift), y1=*(double*)((char*)params+yshift);
|
xue@1
|
657 x2=(x0+xb)/2;
|
xue@1
|
658 ddy2=ddy(x2, params);
|
xue@1
|
659 dy2=*(double*)((char*)params+y1shift), y2=*(double*)((char*)params+yshift);
|
xue@1
|
660 if (y0<=y1 && y0<=y2) xa=x1, ya=y1, dya=dy1, xb=x2, yb=y2, dyb=dy2;
|
xue@1
|
661 else if (y1<=y0 && y1<=y2) {xb=x0, yb=y0, dyb=dy0; x0=x1, y0=y1, dy0=dy1, ddy0=ddy1;}
|
xue@1
|
662 else {xa=x0, ya=y0, dya=dy0; x0=x2, y0=y2, dy0=dy2, ddy0=ddy2;}
|
xue@1
|
663 }
|
xue@1
|
664 while((dya>=0 || dyb<=0) && xb-xa>=epx);
|
xue@1
|
665 }
|
xue@1
|
666 else
|
xue@1
|
667 {
|
xue@1
|
668 //Newton point
|
xue@1
|
669 //dx is the same direction as dy0 if ddy0<0
|
xue@1
|
670 if (fabs(dy0)<epdy) return 0;
|
xue@1
|
671 dx=-dy0/ddy0;
|
xue@1
|
672 x2=x0+dx;
|
xue@1
|
673 if (x2==x0) return 0;
|
xue@1
|
674 //make sure x2 lies in [xa, xb]
|
xue@1
|
675 if (x2<xa) {x2=(xa+x0)/2; dx=x2-x0;}
|
xue@1
|
676 else if (x2>xb) {x2=(xb+x0)/2; dx=x2-x0;}
|
xue@1
|
677 ddy2=ddy(x2, params);
|
xue@1
|
678 y2=*(double*)((char*)params+yshift);
|
xue@1
|
679 //make sure the iteration yields a larger y
|
xue@1
|
680 if (y2>y0)
|
xue@1
|
681 {
|
xue@1
|
682 while (y2>y0)
|
xue@1
|
683 {
|
xue@1
|
684 dx/=2;
|
xue@1
|
685 x2=x0+dx;
|
xue@1
|
686 ddy2=ddy(x2, params);
|
xue@1
|
687 y2=*(double*)((char*)params+yshift);
|
xue@1
|
688 }
|
xue@1
|
689 if (x2==x0) return 0; //dy0<0 but all y(x+(ep)dx)<0
|
xue@1
|
690 }
|
xue@1
|
691 dy2=*(double*)((char*)params+y1shift);
|
xue@1
|
692 if (fabs(dy2)<epdy)
|
xue@1
|
693 {
|
xue@1
|
694 x0=x2;
|
xue@1
|
695 return 0;
|
xue@1
|
696 }
|
xue@1
|
697 else
|
xue@1
|
698 {
|
xue@1
|
699 double fabsdx=fabs(dx);
|
xue@1
|
700 if (fabsdx<epx)
|
xue@1
|
701 {
|
xue@1
|
702 if nsign(dy2, dy0)
|
xue@1
|
703 {
|
xue@1
|
704 x0=x2;
|
xue@1
|
705 return fabsdx;
|
xue@1
|
706 }
|
xue@1
|
707 else
|
xue@1
|
708 {
|
xue@1
|
709 x3=x2+dx;
|
xue@1
|
710 ddy3=ddy(x3, params);
|
xue@1
|
711 dy3=*(double*)((char*)params+y1shift), y3=*(double*)((char*)params+yshift);
|
xue@1
|
712 if (dy3==0) {x0=x3; return 0;}
|
xue@1
|
713 if nsign(dy2, dy3)
|
xue@1
|
714 {
|
xue@1
|
715 x0=x2;
|
xue@1
|
716 return dx;
|
xue@1
|
717 }
|
xue@1
|
718 else if (y3<y2)
|
xue@1
|
719 {
|
xue@1
|
720 x2=x3, y2=y3, dy2=dy3, ddy2=ddy3;
|
xue@1
|
721 }
|
xue@1
|
722 }
|
xue@1
|
723 }
|
xue@1
|
724 }
|
xue@1
|
725 if (x0<x2) xa=x0, ya=y0, dya=dy0;
|
xue@1
|
726 else xb=x0, yb=y0, dyb=dy0;
|
xue@1
|
727
|
xue@1
|
728 x0=x2, y0=y2, dy0=dy2, ddy0=ddy2;
|
xue@1
|
729 }
|
xue@1
|
730
|
xue@1
|
731 iter++;
|
xue@1
|
732 }
|
xue@1
|
733 if (fabs(dy0)<epdy) return 0;
|
xue@1
|
734 else return -1;
|
xue@1
|
735 }//Newton_1d_min
|
xue@1
|
736
|
Chris@5
|
737 /**
|
xue@1
|
738 function Newton1dmax: uses init_Newton_1d_max and Newton_1d_max to search for a local maximum of y.
|
xue@1
|
739
|
xue@1
|
740 By default ddy and dy use the same parameter structure and returns the lower order derivatives at
|
xue@1
|
741 specific offsets within that structure. However, since the same params is used, it's possible the same
|
xue@1
|
742 variable may take different offsets when returned from ddy() and dy().
|
xue@1
|
743
|
xue@1
|
744 On calling xa<x0<xb, and it is assumed that y0>=ya, y0>=yb. If the conditions on y are not met, it
|
xue@1
|
745 returns a negative value showing the order of y0, ya and yb, i.e.
|
xue@1
|
746 -0.5 for ya<=y0<yb
|
xue@1
|
747 -0.6 for yb<=y0<ya
|
xue@1
|
748 -0.7 for y0<ya<yb
|
xue@1
|
749 -0.8 for y0<yb<=ya
|
xue@1
|
750
|
xue@1
|
751 In: function ddy(x) that computes y, y' and y" at x
|
xue@1
|
752 function dy(x) that computes y and y' at x
|
xue@1
|
753 params: additional argument for calling ddy and dy
|
xue@1
|
754 y1shift: the byte shift in params where ddy(x) put the value of y'(x) as double type
|
xue@1
|
755 yshift: the byte shift in params where ddy(x) put the value of y(x) as double type
|
xue@1
|
756 dyshift: the byte shift in params where dy(x) put the value of y(x) as double type
|
xue@1
|
757 x0: inital value of x
|
xue@1
|
758 (xa, xb): initial search range
|
xue@1
|
759 epx, epdy: stopping threshold for x and y', respectively
|
xue@1
|
760 maxiter: maximal number of iterates
|
xue@1
|
761 Out: x0: the maximum
|
xue@1
|
762
|
xue@1
|
763 Returns the error bound for maximum x0, or a negative value if the initial point does not satisfy
|
xue@1
|
764 relevent assumptions.
|
xue@1
|
765 */
|
xue@1
|
766 double Newton1dmax(double& x0, double xa, double xb, double (*ddy)(double, void*), void* params, int y1shift, int yshift, double (*dy)(double, void*), int dyshift, double epx, int maxiter, double epdy, bool checkinput)
|
xue@1
|
767 {
|
xue@1
|
768 double ep=init_Newton_1d_max(x0, xa, xb, dy, params, dyshift, epx);
|
xue@1
|
769 if (ep>epx) ep=Newton_1d_max(x0, xa, xb, ddy, params, y1shift, yshift, epx, maxiter, epdy, true);
|
xue@1
|
770 return ep;
|
xue@1
|
771 }//Newton1dmax
|
xue@1
|
772
|
xue@1
|
773 //function Newton1dmin: same as Newton1dmax, in the other direction.
|
xue@1
|
774 double Newton1dmin(double& x0, double xa, double xb, double (*ddy)(double, void*), void* params, int y1shift, int yshift, double (*dy)(double, void*), int dyshift, double epx, int maxiter, double epdy, bool checkinput)
|
xue@1
|
775 {
|
xue@1
|
776 double ep=init_Newton_1d_min(x0, xa, xb, dy, params, dyshift, epx);
|
xue@1
|
777 if (ep>epx) ep=Newton_1d_min(x0, xa, xb, ddy, params, y1shift, yshift, epx, maxiter, epdy, true);
|
xue@1
|
778 return ep;
|
xue@1
|
779 }//Newton1dmin
|
xue@1
|
780
|
Chris@5
|
781 /**
|
xue@1
|
782 function NewtonMax: looks for the maximum of y within [x0, x1] using Newton method.
|
xue@1
|
783
|
xue@1
|
784 In: function ddy(x) that computes y, y' and y" at x
|
xue@1
|
785 params: additional argument for calling ddy
|
xue@1
|
786 y1shift: the byte shift in params where ddy(x) put the value of y'(x) as double type
|
xue@1
|
787 yshift: the byte shift in params where ddy(x) put the value of y(x) as double type
|
xue@1
|
788 (x0, x1): initial search range
|
xue@1
|
789 epx, epy: stopping threshold for x and y, respectively
|
xue@1
|
790 maxiter: maximal number of iterates
|
xue@1
|
791 Out: maximum x0
|
xue@1
|
792
|
xue@1
|
793 Returns an error bound on maximum x0 if successful, a negative value if not.
|
xue@1
|
794 */
|
xue@1
|
795 double NewtonMax(double& x0, double x1, double (*ddy)(double, void*), void* params, int y1shift, int yshift, double epx, int maxiter, double epy)
|
xue@1
|
796 {
|
xue@1
|
797 double y0, x2, y2, dx, dy0, dy2, dy3, x3, y3;
|
xue@1
|
798 dy0=ddy(x0, params);
|
xue@1
|
799 y0=*(double*)((char*)params+yshift);
|
xue@1
|
800 if (y0==0) return 0;
|
xue@1
|
801
|
xue@1
|
802 int iter=0;
|
xue@1
|
803 while (iter<maxiter)
|
xue@1
|
804 {
|
xue@1
|
805 if (dy0==0) return -1;
|
xue@1
|
806 dx=-y0/dy0;
|
xue@1
|
807
|
xue@1
|
808 x2=x0+dx;
|
xue@1
|
809 dy2=ddy(x2, params);
|
xue@1
|
810 y2=*(double*)((char*)params+yshift);
|
xue@1
|
811 //make sure the iteration yields a y closer to 0
|
xue@1
|
812 while (fabs(y2)>fabs(y0))
|
xue@1
|
813 {
|
xue@1
|
814 dx/=2;
|
xue@1
|
815 x2=x0+dx;
|
xue@1
|
816 dy2=ddy(x2, params);
|
xue@1
|
817 y2=*(double*)((char*)params+yshift);
|
xue@1
|
818 }
|
xue@1
|
819 //the above may fail due to precision problem, i.e. x0+dx=x0
|
xue@1
|
820 if (x2==x0)
|
xue@1
|
821 {
|
xue@1
|
822 if (fabs(y0)<epy) return 0;
|
xue@1
|
823 else
|
xue@1
|
824 {
|
xue@1
|
825 x2=nextdouble(x0, dx);
|
xue@1
|
826 dy2=ddy(x2, params);
|
xue@1
|
827 y2=*(double*)((char*)params+yshift);
|
xue@1
|
828 if (y2==0)
|
xue@1
|
829 {
|
xue@1
|
830 x0=x2;
|
xue@1
|
831 return 0;
|
xue@1
|
832 }
|
xue@1
|
833 else
|
xue@1
|
834 {
|
xue@1
|
835 dx=x2-x0;
|
xue@1
|
836 while (y2==y0)
|
xue@1
|
837 {
|
xue@1
|
838 dx*=2;
|
xue@1
|
839 x2=x0+dx;
|
xue@1
|
840 dy2=ddy(x2, params);
|
xue@1
|
841 y2=*(double*)((char*)params+yshift);
|
xue@1
|
842 }
|
xue@1
|
843 if nsign(y2, y0)
|
xue@1
|
844 return fabs(dx);
|
xue@1
|
845 else if (fabs(y2)>fabs(y0))
|
xue@1
|
846 return -1;
|
xue@1
|
847 }
|
xue@1
|
848 }
|
xue@1
|
849 }
|
xue@1
|
850 else if (y2==0)
|
xue@1
|
851 {
|
xue@1
|
852 x0=x2;
|
xue@1
|
853 return 0;
|
xue@1
|
854 }
|
xue@1
|
855 else
|
xue@1
|
856 {
|
xue@1
|
857 double fabsdx=fabs(dx);
|
xue@1
|
858 if (fabsdx<epx)
|
xue@1
|
859 {
|
xue@1
|
860 if nsign(y2, y0)
|
xue@1
|
861 {
|
xue@1
|
862 x0=x2;
|
xue@1
|
863 return fabsdx;
|
xue@1
|
864 }
|
xue@1
|
865 else
|
xue@1
|
866 {
|
xue@1
|
867 x3=x2+dx;
|
xue@1
|
868 dy3=ddy(x3, params);
|
xue@1
|
869 y3=*(double*)((char*)params+yshift);
|
xue@1
|
870 if (y3==0){x0=x3; return 0;}
|
xue@1
|
871 if nsign(y2, y3)
|
xue@1
|
872 {
|
xue@1
|
873 x0=x2;
|
xue@1
|
874 return dx;
|
xue@1
|
875 }
|
xue@1
|
876 else if (fabs(y3)<fabs(y2))
|
xue@1
|
877 {
|
xue@1
|
878 x2=x3, y2=y3, dy2=dy3;
|
xue@1
|
879 }
|
xue@1
|
880 }
|
xue@1
|
881 }
|
xue@1
|
882 }
|
xue@1
|
883 x0=x2;
|
xue@1
|
884 y0=y2;
|
xue@1
|
885 dy0=dy2;
|
xue@1
|
886 iter++;
|
xue@1
|
887 }
|
xue@1
|
888 if (fabs(y0)<epy) return 0;
|
xue@1
|
889 else return -1;
|
xue@1
|
890 }//NewtonMax
|
xue@1
|
891
|
xue@1
|
892 //---------------------------------------------------------------------------
|
Chris@5
|
893 /**
|
xue@1
|
894 function rootsearchhalf: searches for the root of y(x) in (x1, x2) by half-split method. On calling
|
xue@1
|
895 y(x1) and y(x2) must not have the same sign, or -1 is returned signaling a failure.
|
xue@1
|
896
|
xue@1
|
897 In: function y(x),
|
xue@1
|
898 params: additional arguments for calling y
|
xue@1
|
899 (x1, x2): inital range
|
xue@1
|
900 epx: stopping threshold on x
|
xue@1
|
901 maxiter: maximal number of iterates
|
xue@1
|
902 Out: root x1
|
xue@1
|
903
|
xue@1
|
904 Returns a positive error estimate if a root is found, -1 if not.
|
xue@1
|
905 */
|
xue@1
|
906 double rootsearchhalf(double&x1, double x2, double (*y)(double, void*), void* params, double epx, int maxiter)
|
xue@1
|
907 {
|
xue@1
|
908 if (x1>x2) {double tmp=x1; x1=x2; x2=tmp;}
|
xue@1
|
909 double y1=y(x1, params), y3, x3, ep2=epx*2;
|
xue@1
|
910 if (y1==0) return 0;
|
xue@1
|
911 else
|
xue@1
|
912 {
|
xue@1
|
913 double y2=y(x2, params);
|
xue@1
|
914 if (y2==0)
|
xue@1
|
915 {
|
xue@1
|
916 x1=x2;
|
xue@1
|
917 return 0;
|
xue@1
|
918 }
|
xue@1
|
919 else if nsign(y1, y2) //y1 and y2 have different signs
|
xue@1
|
920 {
|
xue@1
|
921 int iter=0;
|
xue@1
|
922 while (x2-x1>ep2 && iter<maxiter)
|
xue@1
|
923 {
|
xue@1
|
924 x3=(x1+x2)/2;
|
xue@1
|
925 y3=y(x3, params);
|
xue@1
|
926 if (y3==0)
|
xue@1
|
927 {
|
xue@1
|
928 x1=x3;
|
xue@1
|
929 return 0;
|
xue@1
|
930 }
|
xue@1
|
931 else
|
xue@1
|
932 {
|
xue@1
|
933 if nsign(y1, y3) x2=x3, y2=y3;
|
xue@1
|
934 else x1=x3, y1=y3;
|
xue@1
|
935 }
|
xue@1
|
936 iter++;
|
xue@1
|
937 }
|
xue@1
|
938 x1=(x1+x2)/2;
|
xue@1
|
939 return x2-x1;
|
xue@1
|
940 }
|
xue@1
|
941 else
|
xue@1
|
942 return -1;
|
xue@1
|
943 }
|
xue@1
|
944 }//rootsearchhalf
|
xue@1
|
945
|
xue@1
|
946 //---------------------------------------------------------------------------
|
Chris@5
|
947 /**
|
xue@1
|
948 function rootsearchsecant: searches for the root of y(x) in (x1, x2) by secant method.
|
xue@1
|
949
|
xue@1
|
950 In: function y(x),
|
xue@1
|
951 params: additional arguments for calling y
|
xue@1
|
952 (x1, x2): inital range
|
xue@1
|
953 epx: stopping threshold on x
|
xue@1
|
954 maxiter: maximal number of iterates
|
xue@1
|
955 Out: root x1
|
xue@1
|
956
|
xue@1
|
957 Returns a positive error estimate if a root is found, -1 if not.
|
xue@1
|
958 */
|
xue@1
|
959 double rootsearchsecant(double& x1, double x2, double (*y)(double, void*), void* params, double epx, int maxiter)
|
xue@1
|
960 {
|
xue@1
|
961 double y1=y(x1, params), y2=y(x2, params), y3, x3, dx, x0, y0;
|
xue@1
|
962 if (y1==0) return 0;
|
xue@1
|
963 if (y2==0)
|
xue@1
|
964 {
|
xue@1
|
965 x1=x2;
|
xue@1
|
966 return 0;
|
xue@1
|
967 }
|
xue@1
|
968 if (fabs(y1)>fabs(y2))
|
xue@1
|
969 {
|
xue@1
|
970 x3=x1, y3=y1;
|
xue@1
|
971 }
|
xue@1
|
972 else
|
xue@1
|
973 {
|
xue@1
|
974 x3=x2, y3=y2;
|
xue@1
|
975 x2=x1, y2=y1;
|
xue@1
|
976 }
|
xue@1
|
977
|
xue@1
|
978 int iter=0;
|
xue@1
|
979 while (iter<maxiter)
|
xue@1
|
980 {
|
xue@1
|
981 x1=(y3*x2-y2*x3)/(y3-y2);
|
xue@1
|
982 y1=y(x1, params);
|
xue@1
|
983
|
xue@1
|
984 dx=fabs(x1-x2);
|
xue@1
|
985 if (dx<epx)
|
xue@1
|
986 {
|
xue@1
|
987 if nsign(y1, y2)
|
xue@1
|
988 {
|
xue@1
|
989 return dx;
|
xue@1
|
990 }
|
xue@1
|
991 else
|
xue@1
|
992 {
|
xue@1
|
993 x0=2*x1-x2;
|
xue@1
|
994 y0=y(x0, params);
|
xue@1
|
995 if nsign(y0, y1)
|
xue@1
|
996 {
|
xue@1
|
997 return dx;
|
xue@1
|
998 }
|
xue@1
|
999 }
|
xue@1
|
1000 }
|
xue@1
|
1001 x3=x2, y3=y2;
|
xue@1
|
1002 x2=x1, y2=y1;
|
xue@1
|
1003 iter++;
|
xue@1
|
1004 }
|
xue@1
|
1005 return -1;
|
xue@1
|
1006 }//rootsearchsecant
|
xue@1
|
1007
|
Chris@5
|
1008 /**
|
xue@1
|
1009 function rootsearchbisecant: searches for the root of y(x) in (x1, x2) by bi-secant method. On calling
|
xue@1
|
1010 given that y(x1) and y(x2) have different signs. The bound interval (x1, x2) converges to 0, as
|
xue@1
|
1011 compared to standared secant method.
|
xue@1
|
1012
|
xue@1
|
1013 In: function y(x),
|
xue@1
|
1014 params: additional arguments for calling y
|
xue@1
|
1015 (x1, x2): inital range
|
xue@1
|
1016 epx: stopping threshold on x
|
xue@1
|
1017 maxiter: maximal number of iterates
|
xue@1
|
1018 Out: root x1
|
xue@1
|
1019
|
xue@1
|
1020 Returns an estimated error bound if a root is found, a negative value if not.
|
xue@1
|
1021 */
|
xue@1
|
1022 double rootsearchbisecant(double&x1, double x2, double (*y)(double, void*), void* params, double epx, int maxiter)
|
xue@1
|
1023 {
|
xue@1
|
1024 if (x1>x2) {double tmp=x1; x1=x2; x2=tmp;}
|
xue@1
|
1025 double y1=y(x1, params), y3, y4, dy2, dy3, x3, x4, ep2=epx*2;
|
xue@1
|
1026 if (y1==0) return 0;
|
xue@1
|
1027 else
|
xue@1
|
1028 {
|
xue@1
|
1029 double y2=y(x2, params);
|
xue@1
|
1030 if (y2==0)
|
xue@1
|
1031 {
|
xue@1
|
1032 x1=x2;
|
xue@1
|
1033 return 0;
|
xue@1
|
1034 }
|
xue@1
|
1035 else if nsign(y1, y2) //y1 and y2 have different signs
|
xue@1
|
1036 {
|
xue@1
|
1037 int iter=0;
|
xue@1
|
1038 while (x2-x1>ep2 && iter<maxiter)
|
xue@1
|
1039 {
|
xue@1
|
1040 dy2=y1-y2;
|
xue@1
|
1041 x3=(y1*x2-y2*x1)/dy2;
|
xue@1
|
1042 y3=y(x3, params);
|
xue@1
|
1043 if (y3==0)
|
xue@1
|
1044 {
|
xue@1
|
1045 x1=x3;
|
xue@1
|
1046 return 0;
|
xue@1
|
1047 }
|
xue@1
|
1048 else
|
xue@1
|
1049 {
|
xue@1
|
1050 if nsign(y1, y3)
|
xue@1
|
1051 {
|
xue@1
|
1052 dy3=y3-y2;
|
xue@1
|
1053 if (dy3==0 || nsign(dy3, dy2)) x2=x3, y2=y3;
|
xue@1
|
1054 else
|
xue@1
|
1055 {
|
xue@1
|
1056 x4=(y3*x2-y2*x3)/dy3; //x4<x3
|
xue@1
|
1057 if (x4<=x1) x2=x3, y2=y3;
|
xue@1
|
1058 else
|
xue@1
|
1059 {
|
xue@1
|
1060 y4=y(x4, params);
|
xue@1
|
1061 if (y4==0)
|
xue@1
|
1062 {
|
xue@1
|
1063 x1=x4;
|
xue@1
|
1064 return 0;
|
xue@1
|
1065 }
|
xue@1
|
1066 else
|
xue@1
|
1067 {
|
xue@1
|
1068 if nsign(y4, y3) x1=x4, y1=y4, x2=x3, y2=y3;
|
xue@1
|
1069 else x2=x4, y2=y4;
|
xue@1
|
1070 }
|
xue@1
|
1071 }
|
xue@1
|
1072 }
|
xue@1
|
1073 }
|
xue@1
|
1074 else
|
xue@1
|
1075 {
|
xue@1
|
1076 dy3=y1-y3;
|
xue@1
|
1077 if (dy3==0 || nsign(dy3, dy2)) x1=x3, y1=y3;
|
xue@1
|
1078 else
|
xue@1
|
1079 {
|
xue@1
|
1080 x4=(y1*x3-y3*x1)/dy3; //x4>x3
|
xue@1
|
1081 if (x4>=x2) x1=x3, y1=y3;
|
xue@1
|
1082 else
|
xue@1
|
1083 {
|
xue@1
|
1084 y4=y(x4, params);
|
xue@1
|
1085 if (y4==0)
|
xue@1
|
1086 {
|
xue@1
|
1087 x1=x4;
|
xue@1
|
1088 return 0;
|
xue@1
|
1089 }
|
xue@1
|
1090 else
|
xue@1
|
1091 {
|
xue@1
|
1092 if nsign(y4, y3) x1=x3, y1=y3, x2=x4, y2=y4;
|
xue@1
|
1093 else x1=x4, y1=y4;
|
xue@1
|
1094 }
|
xue@1
|
1095 }
|
xue@1
|
1096 }
|
xue@1
|
1097 }
|
xue@1
|
1098 }
|
xue@1
|
1099 iter++;
|
xue@1
|
1100 }
|
xue@1
|
1101 x1=(x1+x2)/2;
|
xue@1
|
1102 return x2-x1;
|
xue@1
|
1103 }
|
xue@1
|
1104 else
|
xue@1
|
1105 return -1;
|
xue@1
|
1106 }
|
xue@1
|
1107 }//rootsearchbisecant
|
xue@1
|
1108
|
xue@1
|
1109 //---------------------------------------------------------------------------
|
Chris@5
|
1110 /**
|
xue@1
|
1111 function Search1Dmax: 1-dimensional maximum search within interval [inf, sup] using 0.618 method
|
xue@1
|
1112
|
xue@1
|
1113 In: function F(x),
|
xue@1
|
1114 params: additional arguments for calling F
|
xue@1
|
1115 (inf, sup): inital range
|
xue@1
|
1116 epx: stopping threshold on x
|
xue@1
|
1117 maxiter: maximal number of iterates
|
xue@1
|
1118 convcheck: specifies if convexity is checked at each interate, false recommended.
|
xue@1
|
1119 Out: x: the maximum
|
xue@1
|
1120 maxresult: the maximal value, if specified
|
xue@1
|
1121
|
xue@1
|
1122 Returns an estimated error bound of x if successful, a negative value if not.
|
xue@1
|
1123 */
|
xue@1
|
1124 double Search1Dmax(double& x, void* params, double (*F)(double, void*), double inf, double sup, double* maxresult, double ep, int maxiter, bool convcheck)
|
xue@1
|
1125 {
|
xue@1
|
1126 double v0618=(sqrt(5.0)-1)/2, iv0618=1-v0618;
|
xue@1
|
1127 double l=sup-inf;
|
xue@1
|
1128 double startsup=sup, startinf=inf, startl=l;
|
xue@1
|
1129 double x1=inf+l*(1-v0618);
|
xue@1
|
1130 double x2=inf+l*v0618;
|
xue@1
|
1131 double result;
|
xue@1
|
1132
|
xue@1
|
1133 double v0=F(inf, params);
|
xue@1
|
1134 double v1=F(x1, params);
|
xue@1
|
1135 double v2=F(x2, params);
|
xue@1
|
1136 double v3=F(sup, params);
|
xue@1
|
1137 int iter=0;
|
xue@1
|
1138
|
xue@1
|
1139 if (v1<v0 && v1<v3 && v2<v0 && v2<v3) goto return1;
|
xue@1
|
1140 if (convcheck && (v1<v0*v0618+v3*iv0618 || v2<v0*iv0618+v3*v0618 || v1<v0*iv0618+v2*v0618 || v2<v1*v0618+v3*iv0618)) goto return1;
|
xue@1
|
1141
|
xue@1
|
1142 while (iter<maxiter)
|
xue@1
|
1143 {
|
xue@1
|
1144 if (l<startl*ep) break;
|
xue@1
|
1145 if (v1>v2)
|
xue@1
|
1146 {
|
xue@1
|
1147 sup=x2, v3=v2;
|
xue@1
|
1148 l=sup-inf;
|
xue@1
|
1149 x2=x1, v2=v1;
|
xue@1
|
1150 x1=inf+l*(1-v0618);
|
xue@1
|
1151 v1=F(x1, params);
|
xue@1
|
1152 if (convcheck && (v1<v0*v0618+v3*iv0618 || v1<v0*iv0618+v2*v0618 || v2<v1*v0618+v3*iv0618)) goto return1;
|
xue@1
|
1153 }
|
xue@1
|
1154 else
|
xue@1
|
1155 {
|
xue@1
|
1156 inf=x1, v0=v1;
|
xue@1
|
1157 l=sup-inf;
|
xue@1
|
1158 x1=x2, v1=v2;
|
xue@1
|
1159 x2=inf+l*v0618;
|
xue@1
|
1160 v2=F(x2, params);
|
xue@1
|
1161 if (convcheck && (v2<v0*iv0618+v3*v0618 || v1<v0*iv0618+v2*v0618 || v2<v1*v0618+v3*iv0618)) goto return1;
|
xue@1
|
1162 }
|
xue@1
|
1163 iter++;
|
xue@1
|
1164 }
|
xue@1
|
1165 if (v1>v2) {x=x1; if (maxresult) maxresult[0]=v1; result=x1-inf;}
|
xue@1
|
1166 else {x=x2; if (maxresult) maxresult[0]=v2; result=sup-x2;}
|
xue@1
|
1167 return result;
|
xue@1
|
1168
|
xue@1
|
1169 return1: //where at some stage v(.) fails to be convex
|
xue@1
|
1170 if (v0>=v1 && v0>=v2 && v0>=v3){x=inf; if (maxresult) maxresult[0]=v0; result=(x==startinf)?(-1):(sup-x1);}
|
xue@1
|
1171 else if (v1>=v0 && v1>=v2 && v1>=v3){x=x1; if (maxresult) maxresult[0]=v1; result=x1-inf;}
|
xue@1
|
1172 else if (v2>=v0 && v2>=v1 && v2>=v3){x=x2; if (maxresult) maxresult[0]=v2; result=sup-x2;}
|
xue@1
|
1173 else {x=sup; if (maxresult) maxresult[0]=v3; result=(x==startsup)?(-1):(x2-inf);}
|
xue@1
|
1174 return result;
|
xue@1
|
1175 }//Search1Dmax
|
xue@1
|
1176
|
Chris@5
|
1177 /**
|
xue@1
|
1178 function Search1DmaxEx: 1-dimensional maximum search using Search1Dmax, but allows setting a $len
|
xue@1
|
1179 argument so that a pre-location of the maximum is performed by sampling (inf, sup) at interval $len,
|
xue@1
|
1180 so that Search1Dmax will work on an interval no longer than 2*len.
|
xue@1
|
1181
|
xue@1
|
1182 In: function F(x),
|
xue@1
|
1183 params: additional arguments for calling F
|
xue@1
|
1184 (inf, sup): inital range
|
xue@1
|
1185 len: rough sampling interval for pre-locating maximum
|
xue@1
|
1186 epx: stopping threshold on x
|
xue@1
|
1187 maxiter: maximal number of iterates
|
xue@1
|
1188 convcheck: specifies if convexity is checked at each interate, false recommended.
|
xue@1
|
1189 Out: x: the maximum
|
xue@1
|
1190 maxresult: the maximal value, if specified
|
xue@1
|
1191
|
xue@1
|
1192 Returns an estimated error bound of x if successful, a negative value if not.
|
xue@1
|
1193 */
|
xue@1
|
1194 double Search1DmaxEx(double& x, void* params, double (*F)(double, void*), double inf, double sup, double* maxresult, double ep, int maxiter, double len)
|
xue@1
|
1195 {
|
xue@1
|
1196 int c=ceil((sup-inf)/len+1);
|
xue@1
|
1197 if (c<8) c=8;
|
xue@1
|
1198 double* vs=new double[c+1], step=(sup-inf)/c;
|
xue@1
|
1199 for (int i=0; i<=c; i++) vs[i]=F(inf+step*i, params);
|
xue@1
|
1200 int max=0; for (int i=1; i<=c; i++) if (vs[max]<vs[i]) max=i;
|
xue@1
|
1201 *maxresult=vs[max]; delete[] vs;
|
xue@1
|
1202 if (max>0 && max<c) return Search1Dmax(x, params, F, inf+step*(max-1), inf+step*(max+1), maxresult, ep, maxiter);
|
xue@1
|
1203 else if (max==0)
|
xue@1
|
1204 {
|
xue@1
|
1205 if (step<=ep) {x=inf; return ep;}
|
xue@1
|
1206 else return Search1DmaxEx(x, params, F, inf, inf+step, maxresult, ep, maxiter, len);
|
xue@1
|
1207 }
|
xue@1
|
1208 else
|
xue@1
|
1209 {
|
xue@1
|
1210 if (step<=ep) {x=sup; return ep;}
|
xue@1
|
1211 else return Search1DmaxEx(x, params, F, sup-step, sup, maxresult, ep, maxiter, len);
|
xue@1
|
1212 }
|
xue@1
|
1213 }//Search1DmaxEx
|
xue@1
|
1214
|
Chris@5
|
1215 /**
|
xue@1
|
1216 function Search1Dmaxfrom: 1-dimensional maximum search from start[] towards start+direct[]. The
|
xue@1
|
1217 function tries to locate an interval ("step") for which the 4-point convexity rule confirms the
|
xue@1
|
1218 existence of a convex maximum.
|
xue@1
|
1219
|
xue@1
|
1220 In: function F(x)
|
xue@1
|
1221 params: additional arguments for calling F
|
xue@1
|
1222 start[dim]: initial x[]
|
xue@1
|
1223 direct[dim]: search direction
|
xue@1
|
1224 step: initial step size
|
xue@1
|
1225 minstep: minimal step size in 4-point scheme between x0 and x3
|
xue@1
|
1226 maxstep: maximal step size in 4-point scheme between x0 and x3
|
xue@1
|
1227 Out: opt[dim]: the maximum, if specified
|
xue@1
|
1228 max: maximal value, if specified
|
xue@1
|
1229
|
xue@1
|
1230 Returns x>0 that locally maximizes F(start+x*direct)
|
xue@1
|
1231 or x|-x>0 that maximizes F(start-x*direct) at calculated points if convex check fails.
|
xue@1
|
1232 */
|
xue@1
|
1233 double Search1Dmaxfrom(void* params, double (*F)(int, double*, void*), int dim, double* start, double* direct,
|
xue@1
|
1234 double step, double minstep, double maxstep, double* opt, double* max, double ep, int maxiter, bool convcheck)
|
xue@1
|
1235 {
|
xue@1
|
1236 bool record=true;
|
xue@1
|
1237 if (!opt)
|
xue@1
|
1238 {
|
xue@1
|
1239 record=false;
|
xue@1
|
1240 opt=new double[dim];
|
xue@1
|
1241 }
|
xue@1
|
1242
|
xue@1
|
1243 double v0618=(sqrt(5.0)-1)/2, iv0618=1-v0618;
|
xue@1
|
1244 double x3=step, x0=0, x1=-1, x2=-1, l=x3-x0;
|
xue@1
|
1245
|
xue@1
|
1246 //prepares x0~x3, v0~v3, x0=0, v2>v3
|
xue@1
|
1247 memcpy(opt, start, sizeof(double)*dim);
|
xue@1
|
1248 double v0=F(dim, opt, params);
|
xue@1
|
1249 MultiAdd(dim, opt, start, direct, x3);
|
xue@1
|
1250 double v1, v2, v3=F(dim, opt, params);
|
xue@1
|
1251 if (v0>=v3)
|
xue@1
|
1252 {
|
xue@1
|
1253 x1=x0+l*iv0618;
|
xue@1
|
1254 MultiAdd(dim, opt, start, direct, x1);
|
xue@1
|
1255 v1=F(dim, opt, params);
|
xue@1
|
1256 x2=x0+l*v0618;
|
xue@1
|
1257 MultiAdd(dim, opt, start, direct, x2);
|
xue@1
|
1258 v2=F(dim, opt, params);
|
xue@1
|
1259 while (l>minstep && v0>v1)
|
xue@1
|
1260 {
|
xue@1
|
1261 x3=x2, v3=v2;
|
xue@1
|
1262 x2=x1, v2=v1;
|
xue@1
|
1263 l=x3-x0;
|
xue@1
|
1264 x1=x0+l*iv0618;
|
xue@1
|
1265 MultiAdd(dim, opt, start, direct, x1);
|
xue@1
|
1266 v1=F(dim, opt, params);
|
xue@1
|
1267 }
|
xue@1
|
1268 if (l<=minstep)
|
xue@1
|
1269 {
|
xue@1
|
1270 if (max) max[0]=v0;
|
xue@1
|
1271 if (!record) delete[] opt;
|
xue@1
|
1272 return x0;
|
xue@1
|
1273 }
|
xue@1
|
1274 }
|
xue@1
|
1275 else//v0<=v3
|
xue@1
|
1276 {
|
xue@1
|
1277 do
|
xue@1
|
1278 {
|
xue@1
|
1279 x1=x2, v1=v2;
|
xue@1
|
1280 x2=x3, v2=v3;
|
xue@1
|
1281 x3+=l*v0618;
|
xue@1
|
1282 l=x3-x0;
|
xue@1
|
1283 MultiAdd(dim, opt, start, direct, x3);
|
xue@1
|
1284 v3=F(dim, opt, params);
|
xue@1
|
1285 } while (v2<v3 && x3-x0<maxstep);
|
xue@1
|
1286 if (l>=maxstep) //increasing all the way to maxstep
|
xue@1
|
1287 {
|
xue@1
|
1288 if (max) max[0]=v3;
|
xue@1
|
1289 if (!record) delete opt;
|
xue@1
|
1290 return x3;
|
xue@1
|
1291 }
|
xue@1
|
1292 else if (x1<0)
|
xue@1
|
1293 {
|
xue@1
|
1294 x1=x0+l*iv0618;
|
xue@1
|
1295 MultiAdd(dim, opt, start, direct, x1);
|
xue@1
|
1296 v1=F(dim, opt, params);
|
xue@1
|
1297 }
|
xue@1
|
1298 }
|
xue@1
|
1299
|
xue@1
|
1300 double oldl=l;
|
xue@1
|
1301
|
xue@1
|
1302 int iter=0;
|
xue@1
|
1303
|
xue@1
|
1304 if (convcheck && (v1<v0*v0618+v3*iv0618 || v2<v0*iv0618+v3*v0618
|
xue@1
|
1305 || v1<v0*iv0618+v2*v0618 || v2<v1*v0618+v3*iv0618)) goto return1;
|
xue@1
|
1306
|
xue@1
|
1307 while (iter<maxiter)
|
xue@1
|
1308 {
|
xue@1
|
1309 if (l<oldl*ep) break;
|
xue@1
|
1310 if (v1>v2)
|
xue@1
|
1311 {
|
xue@1
|
1312 x3=x2, v3=v2;
|
xue@1
|
1313 l=x3-x0;
|
xue@1
|
1314 x2=x1, v2=v1;
|
xue@1
|
1315 x1=x0+l*iv0618;
|
xue@1
|
1316 MultiAdd(dim, opt, start, direct, x1);
|
xue@1
|
1317 v1=F(dim, opt, params);
|
xue@1
|
1318 if (convcheck && (v1<v0*v0618+v3*iv0618 || v1<v0*iv0618+v2*v0618 || v2<v1*v0618+v3*iv0618)) goto return1;
|
xue@1
|
1319 }
|
xue@1
|
1320 else
|
xue@1
|
1321 {
|
xue@1
|
1322 x0=x1, v0=v1;
|
xue@1
|
1323 l=x3-x0;
|
xue@1
|
1324 x1=x2, v1=v2;
|
xue@1
|
1325 x2=x0+l*v0618;
|
xue@1
|
1326 MultiAdd(dim, opt, start, direct, x2);
|
xue@1
|
1327 v2=F(dim, opt, params);
|
xue@1
|
1328 if (convcheck && (v2<v0*iv0618+v3*v0618 || v1<v0*iv0618+v2*v0618 || v2<v1*v0618+v3*iv0618)) goto return1;
|
xue@1
|
1329 }
|
xue@1
|
1330 iter++;
|
xue@1
|
1331 }
|
xue@1
|
1332
|
xue@1
|
1333 if (!record) delete[] opt;
|
xue@1
|
1334
|
xue@1
|
1335 if (max) max[0]=v1;
|
xue@1
|
1336 return x1;
|
xue@1
|
1337
|
xue@1
|
1338 return1:
|
xue@1
|
1339 if (!record) delete[] opt;
|
xue@1
|
1340 if (v1>v0) v0=v1,x0=x1; if (v2>v0) v0=v2, x0=x2; if (v3>v0) v0=v3, x0=x3;
|
xue@1
|
1341 if (max) max[0]=v0;
|
xue@1
|
1342 return -x0;
|
xue@1
|
1343 }//Search1Dmaxfrom
|
xue@1
|
1344
|
xue@1
|
1345
|