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

Last change on this file since 154 was 2, checked in by campagne, 19 years ago

first import
File size: 9.1 KB
Line 
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.