source: GLBFrejus/HEAD/src/glb_min_sup.c@ 151

Last change on this file since 151 was 2, checked in by campagne, 20 years ago

first import
File size: 9.1 KB
RevLine 
[2]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
36static void glb_minimizer_error(char error_text[])
37{
38 glb_fatal(error_text);
39}
40
41double *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
48double **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
64void glb_free_vec(double *v,int nl,int nh)
65{
66 glb_free((char*) (v+nl));
67}
68
69void 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
82typedef struct {
83 int ncom;
84 double *pcom;
85 double *xicom;
86 double (*nrfunc)(double*);
87} glb_min_data;
88
89
90static 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
110static 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
188static 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
250static 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 */
292static double sqrarg;
293#define SQR(a) (sqrarg=(a),sqrarg*sqrarg)
294
295void 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
357void 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
Note: See TracBrowser for help on using the repository browser.