1 | /* GLoBES -- General LOng Baseline Experiment Simulator |
---|
2 | * (C) 2002 - 2004, The GLoBES Team |
---|
3 | * |
---|
4 | * GLoBES is mainly intended for academic purposes. Proper |
---|
5 | * credit must be given if you use GLoBES or parts of it. Please |
---|
6 | * read the section 'Credit' in the README file. |
---|
7 | * |
---|
8 | * This program is free software; you can redistribute it and/or modify |
---|
9 | * it under the terms of the GNU General Public License as published by |
---|
10 | * the Free Software Foundation; either version 2 of the License, or |
---|
11 | * (at your option) any later version. |
---|
12 | * |
---|
13 | * This program is distributed in the hope that it will be useful, |
---|
14 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
15 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
16 | * GNU General Public License for more details. |
---|
17 | * |
---|
18 | * You should have received a copy of the GNU General Public License |
---|
19 | * along with this program; if not, write to the Free Software |
---|
20 | * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
---|
21 | */ |
---|
22 | |
---|
23 | |
---|
24 | |
---|
25 | #include <math.h> |
---|
26 | #include <stdio.h> |
---|
27 | |
---|
28 | /* new */ |
---|
29 | /*#include <stdlib.h> |
---|
30 | #include <math.h> |
---|
31 | #include <string.h> */ |
---|
32 | |
---|
33 | #include "glb_error.h" |
---|
34 | #include "glb_min_sup.h" |
---|
35 | |
---|
36 | static void glb_minimizer_error(char error_text[]) |
---|
37 | { |
---|
38 | glb_fatal(error_text); |
---|
39 | } |
---|
40 | |
---|
41 | double *glb_alloc_vec(int nl,int nh) |
---|
42 | { |
---|
43 | double *v; |
---|
44 | v=(double *)glb_malloc((unsigned) (nh-nl+1)*sizeof(double)); |
---|
45 | return v-nl; |
---|
46 | } |
---|
47 | |
---|
48 | double **glb_alloc_mat(int nrl,int nrh, int ncl,int nch) |
---|
49 | { |
---|
50 | int i; |
---|
51 | double **m; |
---|
52 | |
---|
53 | m=(double **) glb_malloc((unsigned) (nrh-nrl+1)*sizeof(double*)); |
---|
54 | m -= nrl; |
---|
55 | for(i=nrl;i<=nrh;i++) |
---|
56 | { |
---|
57 | m[i]=(double *) glb_malloc((unsigned) |
---|
58 | (nch-ncl+1)*sizeof(double)); |
---|
59 | m[i] -= ncl; |
---|
60 | } |
---|
61 | return m; |
---|
62 | } |
---|
63 | |
---|
64 | void glb_free_vec(double *v,int nl,int nh) |
---|
65 | { |
---|
66 | glb_free((char*) (v+nl)); |
---|
67 | } |
---|
68 | |
---|
69 | void glb_free_mat(double **m,int nrl,int nrh,int ncl,int nch) |
---|
70 | { |
---|
71 | int i; |
---|
72 | for(i=nrh;i>=nrl;i--) glb_free((char*) (m[i]+ncl)); |
---|
73 | glb_free((char*) (m+nrl)); |
---|
74 | } |
---|
75 | |
---|
76 | |
---|
77 | |
---|
78 | |
---|
79 | #define TOL 2.0e-4 |
---|
80 | |
---|
81 | |
---|
82 | typedef struct { |
---|
83 | int ncom; |
---|
84 | double *pcom; |
---|
85 | double *xicom; |
---|
86 | double (*nrfunc)(double*); |
---|
87 | } glb_min_data; |
---|
88 | |
---|
89 | |
---|
90 | static double one_dim_projection(double x,glb_min_data *in) |
---|
91 | { |
---|
92 | int j; |
---|
93 | double f,*xt; |
---|
94 | |
---|
95 | xt=glb_alloc_vec(1,in->ncom); |
---|
96 | for (j=1;j<=in->ncom;j++) xt[j]=in->pcom[j]+x*in->xicom[j]; |
---|
97 | f=in->nrfunc(xt); |
---|
98 | glb_free_vec(xt,1,in->ncom); |
---|
99 | return f; |
---|
100 | } |
---|
101 | |
---|
102 | |
---|
103 | |
---|
104 | #define ITMAX 100 |
---|
105 | #define CGOLD 0.3819660 |
---|
106 | #define ZEPS 1.0e-10 |
---|
107 | #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a)) |
---|
108 | #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); |
---|
109 | |
---|
110 | static double glb_brent_min(double ax,double bx,double cx, |
---|
111 | double (*f)(double,glb_min_data*), |
---|
112 | double tol,double *xmin,glb_min_data *in) |
---|
113 | { |
---|
114 | int iter; |
---|
115 | double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm; |
---|
116 | double e=0.0; |
---|
117 | |
---|
118 | a=((ax < cx) ? ax : cx); |
---|
119 | b=((ax > cx) ? ax : cx); |
---|
120 | x=w=v=bx; |
---|
121 | fw=fv=fx=(*f)(x,in); |
---|
122 | for (iter=1;iter<=ITMAX;iter++) { |
---|
123 | xm=0.5*(a+b); |
---|
124 | tol2=2.0*(tol1=tol*fabs(x)+ZEPS); |
---|
125 | if (fabs(x-xm) <= (tol2-0.5*(b-a))) { |
---|
126 | *xmin=x; |
---|
127 | return fx; |
---|
128 | } |
---|
129 | if (fabs(e) > tol1) { |
---|
130 | r=(x-w)*(fx-fv); |
---|
131 | q=(x-v)*(fx-fw); |
---|
132 | p=(x-v)*q-(x-w)*r; |
---|
133 | q=2.0*(q-r); |
---|
134 | if (q > 0.0) p = -p; |
---|
135 | q=fabs(q); |
---|
136 | etemp=e; |
---|
137 | e=d; |
---|
138 | if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x)) |
---|
139 | d=CGOLD*(e=(x >= xm ? a-x : b-x)); |
---|
140 | else { |
---|
141 | d=p/q; |
---|
142 | u=x+d; |
---|
143 | if (u-a < tol2 || b-u < tol2) |
---|
144 | d=SIGN(tol1,xm-x); |
---|
145 | } |
---|
146 | } else { |
---|
147 | d=CGOLD*(e=(x >= xm ? a-x : b-x)); |
---|
148 | } |
---|
149 | u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d)); |
---|
150 | fu=(*f)(u,in); |
---|
151 | if (fu <= fx) { |
---|
152 | if (u >= x) a=x; else b=x; |
---|
153 | SHFT(v,w,x,u) |
---|
154 | SHFT(fv,fw,fx,fu) |
---|
155 | } else { |
---|
156 | if (u < x) a=u; else b=u; |
---|
157 | if (fu <= fw || w == x) { |
---|
158 | v=w; |
---|
159 | w=u; |
---|
160 | fv=fw; |
---|
161 | fw=fu; |
---|
162 | } else if (fu <= fv || v == x || v == w) { |
---|
163 | v=u; |
---|
164 | fv=fu; |
---|
165 | } |
---|
166 | } |
---|
167 | } |
---|
168 | glb_minimizer_error("glb_brent_min"); |
---|
169 | *xmin=x; |
---|
170 | return fx; |
---|
171 | } |
---|
172 | |
---|
173 | #undef ITMAX |
---|
174 | #undef CGOLD |
---|
175 | #undef ZEPS |
---|
176 | #undef SIGN |
---|
177 | |
---|
178 | |
---|
179 | |
---|
180 | |
---|
181 | #define GOLD 1.618034 |
---|
182 | #define GLIMIT 100.0 |
---|
183 | #define TINY 1.0e-20 |
---|
184 | #define MAX(a,b) ((a) > (b) ? (a) : (b)) |
---|
185 | #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a)) |
---|
186 | #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); |
---|
187 | |
---|
188 | static void bracket(double *ax,double *bx,double *cx,double *fa,double *fb, |
---|
189 | double *fc, double (*func)(double,glb_min_data*),glb_min_data *in) |
---|
190 | { |
---|
191 | double ulim,u,r,q,fu,dum; |
---|
192 | |
---|
193 | *fa=(*func)(*ax,in); |
---|
194 | *fb=(*func)(*bx,in); |
---|
195 | if (*fb > *fa) { |
---|
196 | SHFT(dum,*ax,*bx,dum) |
---|
197 | SHFT(dum,*fb,*fa,dum) |
---|
198 | } |
---|
199 | *cx=(*bx)+GOLD*(*bx-*ax); |
---|
200 | *fc=(*func)(*cx,in); |
---|
201 | while (*fb > *fc) { |
---|
202 | r=(*bx-*ax)*(*fb-*fc); |
---|
203 | q=(*bx-*cx)*(*fb-*fa); |
---|
204 | u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/ |
---|
205 | (2.0*SIGN(MAX(fabs(q-r),TINY),q-r)); |
---|
206 | ulim=(*bx)+GLIMIT*(*cx-*bx); |
---|
207 | if ((*bx-u)*(u-*cx) > 0.0) { |
---|
208 | fu=(*func)(u,in); |
---|
209 | if (fu < *fc) { |
---|
210 | *ax=(*bx); |
---|
211 | *bx=u; |
---|
212 | *fa=(*fb); |
---|
213 | *fb=fu; |
---|
214 | return; |
---|
215 | } else if (fu > *fb) { |
---|
216 | *cx=u; |
---|
217 | *fc=fu; |
---|
218 | return; |
---|
219 | } |
---|
220 | u=(*cx)+GOLD*(*cx-*bx); |
---|
221 | fu=(*func)(u,in); |
---|
222 | } else if ((*cx-u)*(u-ulim) > 0.0) { |
---|
223 | fu=(*func)(u,in); |
---|
224 | if (fu < *fc) { |
---|
225 | SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx)) |
---|
226 | SHFT(*fb,*fc,fu,(*func)(u,in)) |
---|
227 | } |
---|
228 | } else if ((u-ulim)*(ulim-*cx) >= 0.0) { |
---|
229 | u=ulim; |
---|
230 | fu=(*func)(u,in); |
---|
231 | } else { |
---|
232 | u=(*cx)+GOLD*(*cx-*bx); |
---|
233 | fu=(*func)(u,in); |
---|
234 | } |
---|
235 | SHFT(*ax,*bx,*cx,u) |
---|
236 | SHFT(*fa,*fb,*fc,fu) |
---|
237 | } |
---|
238 | } |
---|
239 | |
---|
240 | |
---|
241 | |
---|
242 | #undef GOLD |
---|
243 | #undef GLIMIT |
---|
244 | #undef TINY |
---|
245 | #undef MAX |
---|
246 | #undef SIGN |
---|
247 | #undef SHFT |
---|
248 | |
---|
249 | |
---|
250 | static void line_minimization(double p[],double xi[],int n,double *fret, |
---|
251 | double (*func)(double*), |
---|
252 | glb_min_data *in) |
---|
253 | { |
---|
254 | int j; |
---|
255 | double xx,xmin,fx,fb,fa,bx,ax; |
---|
256 | int ncom=0; |
---|
257 | double *pcom=0,*xicom=0,(*nrfunc)(double*); |
---|
258 | |
---|
259 | ncom=n; |
---|
260 | pcom=glb_alloc_vec(1,n); |
---|
261 | xicom=glb_alloc_vec(1,n); |
---|
262 | nrfunc=func; |
---|
263 | for (j=1;j<=n;j++) { |
---|
264 | pcom[j]=p[j]; |
---|
265 | xicom[j]=xi[j]; |
---|
266 | } |
---|
267 | |
---|
268 | in->ncom=n; |
---|
269 | in->pcom=pcom; |
---|
270 | in->xicom=xicom; |
---|
271 | in->nrfunc=func; |
---|
272 | ax=0.0; |
---|
273 | xx=1.0; |
---|
274 | bx=2.0; |
---|
275 | bracket(&ax,&xx,&bx,&fa,&fx,&fb,one_dim_projection,in); |
---|
276 | *fret=glb_brent_min(ax,xx,bx,one_dim_projection,TOL,&xmin,in); |
---|
277 | for (j=1;j<=n;j++) { |
---|
278 | xi[j] *= xmin; |
---|
279 | p[j] += xi[j]; |
---|
280 | } |
---|
281 | glb_free_vec(xicom,1,n); |
---|
282 | glb_free_vec(pcom,1,n); |
---|
283 | } |
---|
284 | |
---|
285 | |
---|
286 | |
---|
287 | |
---|
288 | |
---|
289 | |
---|
290 | |
---|
291 | #define ITMAX 10000 /* was 500 */ |
---|
292 | static double sqrarg; |
---|
293 | #define SQR(a) (sqrarg=(a),sqrarg*sqrarg) |
---|
294 | |
---|
295 | void glb_powell(double p[],double **xi,int n, |
---|
296 | double ftol,int *iter,double *fret, |
---|
297 | double (*func)(double*)) |
---|
298 | { |
---|
299 | int i,ibig,j; |
---|
300 | double t,fptt,fp,del; |
---|
301 | double *pt,*ptt,*xit; |
---|
302 | glb_min_data in; |
---|
303 | pt=glb_alloc_vec(1,n); |
---|
304 | ptt=glb_alloc_vec(1,n); |
---|
305 | xit=glb_alloc_vec(1,n); |
---|
306 | *fret=(*func)(p); |
---|
307 | for (j=1;j<=n;j++) pt[j]=p[j]; |
---|
308 | for (*iter=1;;(*iter)++) { |
---|
309 | fp=(*fret); |
---|
310 | ibig=0; |
---|
311 | del=0.0; |
---|
312 | for (i=1;i<=n;i++) { |
---|
313 | for (j=1;j<=n;j++) xit[j]=xi[j][i]; |
---|
314 | fptt=(*fret); |
---|
315 | line_minimization(p,xit,n,fret,func,&in); |
---|
316 | if (fabs(fptt-(*fret)) > del) { |
---|
317 | del=fabs(fptt-(*fret)); |
---|
318 | ibig=i; |
---|
319 | } |
---|
320 | } |
---|
321 | if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { |
---|
322 | glb_free_vec(xit,1,n); |
---|
323 | glb_free_vec(ptt,1,n); |
---|
324 | glb_free_vec(pt,1,n); |
---|
325 | return; |
---|
326 | } |
---|
327 | if (*iter == ITMAX) |
---|
328 | { /* modified by WW */ |
---|
329 | /* glb_minimizer_error("empty (1)"); */ |
---|
330 | glb_free_vec(xit,1,n); |
---|
331 | glb_free_vec(ptt,1,n); |
---|
332 | glb_free_vec(pt,1,n); |
---|
333 | (*fret) = 999999.0; |
---|
334 | printf("At Point: "); |
---|
335 | for (j=1;j<=n;j++) printf("%g ",p[j]); |
---|
336 | printf("\n"); |
---|
337 | glb_error("Minimizer error: no convergence (1)"); |
---|
338 | return; |
---|
339 | } |
---|
340 | for (j=1;j<=n;j++) { |
---|
341 | ptt[j]=2.0*p[j]-pt[j]; |
---|
342 | xit[j]=p[j]-pt[j]; |
---|
343 | pt[j]=p[j]; |
---|
344 | } |
---|
345 | fptt=(*func)(ptt); |
---|
346 | if (fptt < fp) { |
---|
347 | t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); |
---|
348 | if (t < 0.0) { |
---|
349 | line_minimization(p,xit,n,fret,func,&in); |
---|
350 | for (j=1;j<=n;j++) xi[j][ibig]=xit[j]; |
---|
351 | } |
---|
352 | } |
---|
353 | } |
---|
354 | } |
---|
355 | |
---|
356 | |
---|
357 | void glb_powell2(double p[],double **xi,int n,double ftol, |
---|
358 | int *iter,double *fret, |
---|
359 | double (*func)(double*)) |
---|
360 | { |
---|
361 | int i,ibig,j; |
---|
362 | double t,fptt,fp,del; |
---|
363 | double *pt,*ptt,*xit; |
---|
364 | glb_min_data in; |
---|
365 | pt=glb_alloc_vec(1,n); |
---|
366 | ptt=glb_alloc_vec(1,n); |
---|
367 | xit=glb_alloc_vec(1,n); |
---|
368 | *fret=(*func)(p); |
---|
369 | for (j=1;j<=n;j++) pt[j]=p[j]; |
---|
370 | for (*iter=1;;(*iter)++) { |
---|
371 | fp=(*fret); |
---|
372 | ibig=0; |
---|
373 | del=0.0; |
---|
374 | for (i=1;i<=n;i++) { |
---|
375 | for (j=1;j<=n;j++) xit[j]=xi[j][i]; |
---|
376 | fptt=(*fret); |
---|
377 | line_minimization(p,xit,n,fret,func,&in); |
---|
378 | if (fabs(fptt-(*fret)) > del) { |
---|
379 | del=fabs(fptt-(*fret)); |
---|
380 | ibig=i; |
---|
381 | } |
---|
382 | } |
---|
383 | if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { |
---|
384 | glb_free_vec(xit,1,n); |
---|
385 | glb_free_vec(ptt,1,n); |
---|
386 | glb_free_vec(pt,1,n); |
---|
387 | return; |
---|
388 | } |
---|
389 | if (*iter == ITMAX) |
---|
390 | { /* modified by WW */ |
---|
391 | /* glb_minimizer_error("empty (2)"); */ |
---|
392 | glb_free_vec(xit,1,n); |
---|
393 | glb_free_vec(ptt,1,n); |
---|
394 | glb_free_vec(pt,1,n); |
---|
395 | (*fret) = 999999.0; |
---|
396 | printf("At Point: "); |
---|
397 | for (j=1;j<=n;j++) printf("%g ",p[j]); |
---|
398 | printf("\n"); |
---|
399 | glb_error("Minimizer error: no convergence (2)"); |
---|
400 | return; |
---|
401 | } |
---|
402 | |
---|
403 | for (j=1;j<=n;j++) { |
---|
404 | ptt[j]=2.0*p[j]-pt[j]; |
---|
405 | xit[j]=p[j]-pt[j]; |
---|
406 | pt[j]=p[j]; |
---|
407 | } |
---|
408 | fptt=(*func)(ptt); |
---|
409 | if (fptt < fp) { |
---|
410 | t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); |
---|
411 | if (t < 0.0) { |
---|
412 | line_minimization(p,xit,n,fret,func,&in); |
---|
413 | for (j=1;j<=n;j++) xi[j][ibig]=xit[j]; |
---|
414 | } |
---|
415 | } |
---|
416 | } |
---|
417 | } |
---|
418 | |
---|
419 | |
---|
420 | |
---|
421 | #undef ITMAX |
---|
422 | #undef SQR |
---|
423 | #undef TOL |
---|