1 //
2 // C++ Implementation: optimization
3 //
4 // Description:
5 //
6 //
7 // Author: BUI Quang Minh, Steffen Klaere, Arndt von Haeseler <minh.bui@univie.ac.at>, (C) 2008
8 //
9 // Copyright: See COPYING file that comes with this distribution
10 //
11 //
12 #include "optimization.h"
13 #include <cmath>
14 #include <stdio.h>
15 #include <stdlib.h>
16 #include <iostream>
17 #include "lbfgsb/lbfgsb_new.h"
18 #include "tools.h"
19
20
21 using namespace std;
22
23 const double ERROR_X = 1.0e-4;
24
25 double ran1(long *idum);
26 double *new_vector(long nl, long nh);
27 void free_vector(double *v, long nl, long nh);
28 double **new_matrix(long nrl, long nrh, long ncl, long nch);
29 void free_matrix(double **m, long nrl, long nrh, long ncl, long nch);
30 void fixBound(double x[], double lower[], double upper[], int n);
31
32 #define NR_END 1
33 #define FREE_ARG char*
34
35 #define GET_PSUM \
36 for (n=1;n<=ndim;n++) {\
37 for (sum=0.0,m=1;m<=mpts;m++) sum += p[m][n];\
38 psum[n]=sum;}
39
40
41 /*
42 #define IA 16807
43 #define IM 2147483647
44 #define AM (1.0/IM)
45 #define IQ 127773
46 #define IR 2836
47 #define NTAB 32
48 #define NDIV (1+(IM-1)/NTAB)
49 #define EPS 1.2e-7
50 #define RNMX (1.0-EPS)
51
52 double ran1(long *idum) {
53 int j;
54 long k;
55 static long iy=0;
56 static long iv[NTAB];
57 double temp;
58
59 if (*idum <= 0 || !iy) {
60 if (-(*idum) < 1) *idum=1;
61 else *idum = -(*idum);
62 for (j=NTAB+7;j>=0;j--) {
63 k=(*idum)/IQ;
64 *idum=IA*(*idum-k*IQ)-IR*k;
65 if (*idum < 0) *idum += IM;
66 if (j < NTAB) iv[j] = *idum;
67 }
68 iy=iv[0];
69 }
70 k=(*idum)/IQ;
71 *idum=IA*(*idum-k*IQ)-IR*k;
72 if (*idum < 0) *idum += IM;
73 j=iy/NDIV;
74 iy=iv[j];
75 iv[j] = *idum;
76 if ((temp=AM*iy) > RNMX) return RNMX;
77 else return temp;
78 }
79 #undef IA
80 #undef IM
81 #undef AM
82 #undef IQ
83 #undef IR
84 #undef NTAB
85 #undef NDIV
86 #undef EPS
87 #undef RNMX
88 */
89
90 long idum = 123456;
91 double tt;
92
nrerror(const char * error_text)93 void nrerror(const char *error_text)
94 /* Numerical Recipes standard error handler */
95 {
96 cerr << "NUMERICAL ERROR: " << error_text << endl;
97 //exit(1);
98 abort();
99 }
100
new_vector(long nl,long nh)101 double *new_vector(long nl, long nh)
102 /* allocate a double vector with subscript range v[nl..nh] */
103 {
104 double *v;
105
106 v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
107 if (!v) nrerror("allocation failure in vector()");
108 return v-nl+NR_END;
109 }
110
new_matrix(long nrl,long nrh,long ncl,long nch)111 double **new_matrix(long nrl, long nrh, long ncl, long nch)
112 /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
113 {
114 long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
115 double **m;
116
117 /* allocate pointers to rows */
118 m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
119 if (!m) nrerror("allocation failure 1 in matrix()");
120 m += NR_END;
121 m -= nrl;
122
123 /* allocate rows and set pointers to them */
124 m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
125 if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
126 m[nrl] += NR_END;
127 m[nrl] -= ncl;
128
129 for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
130
131 /* return pointer to array of pointers to rows */
132 return m;
133 }
134
135
free_vector(double * v,long nl,long nh)136 void free_vector(double *v, long nl, long nh)
137 /* free a double vector allocated with vector() */
138 {
139 free((FREE_ARG) (v+nl-NR_END));
140 }
141
free_matrix(double ** m,long nrl,long nrh,long ncl,long nch)142 void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
143 /* free a double matrix allocated by dmatrix() */
144 {
145 free((FREE_ARG) (m[nrl]+ncl-NR_END));
146 free((FREE_ARG) (m+nrl-NR_END));
147 }
148
fixBound(double x[],double lower[],double upper[],int n)149 void fixBound(double x[], double lower[], double upper[], int n) {
150 for (int i = 1; i <= n; i++) {
151 if (x[i] < lower[i])
152 x[i] = lower[i];
153 else if (x[i] > upper[i])
154 x[i] = upper[i];
155 }
156 }
157
158 /**********************************************
159 Optimization routines
160 **********************************************/
Optimization()161 Optimization::Optimization()
162 {
163 }
164
165
~Optimization()166 Optimization::~Optimization()
167 {
168 }
169
170 /*****************************************************
171 One dimensional optimization with Brent method
172 *****************************************************/
173
174 #define ITMAX 100
175 #define CGOLD 0.3819660
176 #define GOLD 1.618034
177 #define GLIMIT 100.0
178 #define TINY 1.0e-20
179 #define ZEPS 1.0e-10
180 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
181 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
182
183 /* Brents method in one dimension */
brent_opt(double ax,double bx,double cx,double tol,double * foptx,double * f2optx,double fax,double fbx,double fcx)184 double Optimization::brent_opt (double ax, double bx, double cx, double tol,
185 double *foptx, double *f2optx, double fax, double fbx, double fcx) {
186 int iter;
187 double a,b,d=0,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
188 double xw,wv,vx;
189 double e=0.0;
190
191 a=(ax < cx ? ax : cx);
192 b=(ax > cx ? ax : cx);
193 x=bx;
194 fx=fbx;
195 if (fax < fcx) {
196 w=ax;
197 fw=fax;
198 v=cx;
199 fv=fcx;
200 } else {
201 w=cx;
202 fw=fcx;
203 v=ax;
204 fv=fax;
205 }
206
207 for (iter=1;iter<=ITMAX;iter++) {
208 xm=0.5*(a+b);
209 tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
210 if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
211 *foptx = fx;
212 xw = x-w;
213 wv = w-v;
214 vx = v-x;
215 *f2optx = 2.0*(fv*xw + fx*wv + fw*vx)/
216 (v*v*xw + x*x*wv + w*w*vx);
217 return x;
218 }
219
220 if (fabs(e) > tol1) {
221 r=(x-w)*(fx-fv);
222 q=(x-v)*(fx-fw);
223 p=(x-v)*q-(x-w)*r;
224 q=2.0*(q-r);
225 if (q > 0.0)
226 p = -p;
227 q=fabs(q);
228 etemp=e;
229 e=d;
230 if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
231 d=CGOLD*(e=(x >= xm ? a-x : b-x));
232 else {
233 d=p/q;
234 u=x+d;
235 if (u-a < tol2 || b-u < tol2)
236 d=SIGN(tol1,xm-x);
237 }
238 } else {
239 d=CGOLD*(e=(x >= xm ? a-x : b-x));
240 }
241
242 u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
243 fu=computeFunction(u);
244 if (fu <= fx) {
245 if (u >= x)
246 a=x;
247 else
248 b=x;
249
250 SHFT(v,w,x,u)
251 SHFT(fv,fw,fx,fu)
252 } else {
253 if (u < x)
254 a=u;
255 else
256 b=u;
257 if (fu <= fw || w == x) {
258 v=w;
259 w=u;
260 fv=fw;
261 fw=fu;
262 } else
263 if (fu <= fv || v == x || v == w) {
264 v=u;
265 fv=fu;
266 }
267 }
268 }
269
270 *foptx = fx;
271 xw = x-w;
272 wv = w-v;
273 vx = v-x;
274 *f2optx = 2.0*(fv*xw + fx*wv + fw*vx)/(v*v*xw + x*x*wv + w*w*vx);
275
276 return x;
277 }
278
279 #undef ITMAX
280 #undef CGOLD
281 #undef ZEPS
282 #undef SHFT
283 #undef SIGN
284 #undef GOLD
285 #undef GLIMIT
286 #undef TINY
287
288
minimizeOneDimen(double xmin,double xguess,double xmax,double tolerance,double * fx,double * f2x)289 double Optimization::minimizeOneDimen(double xmin, double xguess, double xmax, double tolerance, double *fx, double *f2x) {
290 double eps, optx, ax, bx, cx, fa, fb, fc;
291 //int converged; /* not converged error flag */
292
293 /* first attempt to bracketize minimum */
294 if (xguess < xmin) xguess = xmin;
295 if (xguess > xmax) xguess = xmax;
296 eps = xguess*tolerance*50.0;
297 ax = xguess - eps;
298 if (ax < xmin) ax = xmin;
299 bx = xguess;
300 cx = xguess + eps;
301 if (cx > xmax) cx = xmax;
302
303 /* check if this works */
304 // compute fb first to save some computation, if any
305 fb = computeFunction(bx);
306 fa = computeFunction(ax);
307 fc = computeFunction(cx);
308
309 /* if it works use these borders else be conservative */
310 if ((fa < fb) || (fc < fb)) {
311 if (ax != xmin) fa = computeFunction(xmin);
312 if (cx != xmax) fc = computeFunction(xmax);
313 ax = xmin;
314 cx = xmax;
315 }
316 /*
317 const int MAX_ROUND = 10;
318 for (i = 0; ((fa < fb-tolerance) || (fc < fb-tolerance)) && (i<MAX_ROUND); i++) {
319 // brent method require that fb is smaller than both fa and fc
320 // find some random values until fb achieve this
321 bx = (((double)rand()) / RAND_MAX)*(cx-ax) + ax;
322 fb = computeFunction(bx);
323 }*/
324
325 /*
326 if ((fa < fb) || (fc < fb)) {
327 if (fa < fc) { bx = ax; fb = fa; } else { bx = cx; fb = fc; }
328 //cout << "WARNING: Initial value for Brent method is set at bound " << bx << endl;
329 }*/
330 // optx = brent_opt(xmin, xguess, xmax, tolerance, fx, f2x, fa, fb, fc);
331 //} else
332 optx = brent_opt(ax, bx, cx, tolerance, fx, f2x, fa, fb, fc);
333 if (*fx > fb) // if worse, return initial value
334 {
335 *fx = computeFunction(bx);
336 return bx;
337 }
338
339 return optx; /* return optimal x */
340 }
341
minimizeOneDimenSafeMode(double xmin,double xguess,double xmax,double tolerance,double * f)342 double Optimization::minimizeOneDimenSafeMode(double xmin, double xguess, double xmax, double tolerance, double *f)
343 {
344 double ferror;
345 double optx = minimizeOneDimen(xmin, xguess, xmax, tolerance, f, &ferror);
346 double fnew;
347 // check value at the boundary
348 if ((optx < xmax) && (fnew = computeFunction(xmax)) <= *f+tolerance) {
349 //if (verbose_mode >= VB_MAX)
350 //cout << "Note from Newton safe mode: " << optx << " (" << f << ") -> " << xmax << " ("<< fnew << ")" << endl;
351 optx = xmax;
352 *f = fnew;
353 }
354 if ((optx > xmin) && (fnew = computeFunction(xmin)) <= *f+tolerance) {
355 //if (verbose_mode >= VB_MAX)
356 //cout << "Note from Newton safe mode: " << optx << " -> " << xmin << endl;
357 optx = xmin;
358 *f = fnew;
359 }
360 return optx;
361 }
362
363 /*****************************************************
364 One dimensional optimization with Newton Raphson
365 only applicable if 1st and 2nd derivatives are easy to compute
366 *****************************************************/
367
368
minimizeNewtonSafeMode(double xmin,double xguess,double xmax,double tolerance,double & f)369 double Optimization::minimizeNewtonSafeMode(double xmin, double xguess, double xmax, double tolerance, double &f)
370 {
371 double optx = minimizeNewton(xmin, xguess, xmax, tolerance, f);
372 double fnew;
373 // check value at the boundary
374 if ((optx < xmax) && (fnew = computeFunction(xmax)) <= f+tolerance) {
375 //if (verbose_mode >= VB_MAX)
376 //cout << "Note from Newton safe mode: " << optx << " (" << f << ") -> " << xmax << " ("<< fnew << ")" << endl;
377 optx = xmax;
378 f = fnew;
379 }
380 if ((optx > xmin) && (fnew = computeFunction(xmin)) <= f+tolerance) {
381 //if (verbose_mode >= VB_MAX)
382 //cout << "Note from Newton safe mode: " << optx << " -> " << xmin << endl;
383 optx = xmin;
384 f = fnew;
385 }
386 return optx;
387 }
388
minimizeNewton(double x1,double xguess,double x2,double xacc,double & d2l,int maxNRStep)389 double Optimization::minimizeNewton(double x1, double xguess, double x2, double xacc, double &d2l, int maxNRStep)
390 {
391 int j;
392 double df,dx,dxold,f;
393 double temp,xh,xl,rts, rts_old, xinit;
394
395 rts = xguess;
396 if (rts < x1) rts = x1;
397 if (rts > x2) rts = x2;
398 xinit = xguess;
399 // finit = fold = fm = computeFuncDerv(rts,f,df);
400 computeFuncDerv(rts,f,df);
401 d2l = df;
402 if (!isfinite(f) || !isfinite(df)) {
403 nrerror("Wrong computeFuncDerv");
404 }
405 if (df >= 0.0 && fabs(f) < xacc) return rts;
406 if (f < 0.0) {
407 xl = rts;
408 xh = x2;
409 } else {
410 xh = rts;
411 xl = x1;
412 }
413
414 dx=dxold=fabs(xh-xl);
415 for (j=1;j<=maxNRStep;j++) {
416 rts_old = rts;
417 if (
418 (df <= 0.0) // function is concave
419 // || (fm > fold + xacc) // increasing
420 || (((rts-xh)*df-f)*((rts-xl)*df-f) >= 0.0) // out of bound
421 //|| (fabs(2.0*f) > fabs(dxold*df)) // converge too slow
422 ) {
423 dxold=dx;
424 dx=0.5*(xh-xl);
425 rts=xl+dx;
426 d2l = df;
427 if (xl == rts) return rts;
428 } else {
429 dxold=dx;
430 dx=f/df;
431 temp=rts;
432 rts -= dx;
433 d2l = df;
434 if (temp == rts) return rts;
435 }
436 if (fabs(dx) < xacc || (j == maxNRStep)) {
437 // if (fm > finit) {
438 // // happen in rare cases that it is worse than starting point: revert init value
439 // fm = computeFunction(xinit);
440 // return xinit;
441 // }
442 return rts_old;
443 // return rts;
444 }
445 // fold = fm;
446 // fm = computeFuncDerv(rts,f,df);
447 computeFuncDerv(rts,f,df);
448 if (!isfinite(f) || !isfinite(df)) nrerror("Wrong computeFuncDerv");
449 if (df > 0.0 && fabs(f) < xacc) {
450 d2l = df;
451 // if (fm > finit) {
452 // // happen in rare cases that it is worse than starting point: revert init value
453 // fm = computeFunction(xinit);
454 // return xinit;
455 // }
456 return rts;
457 }
458 if (f < 0.0)
459 xl=rts;
460 else if (f > 0.0)
461 xh=rts;
462 }
463 nrerror("Maximum number of iterations exceeded in minimizeNewton");
464 d2l = 0.0;
465 return 0.0;
466 }
467
minimizeNewton(double x1,double xguess,double x2,double xacc,int maxNRStep)468 double Optimization::minimizeNewton(double x1, double xguess, double x2, double xacc, int maxNRStep)
469 {
470 double var;
471 double optx = minimizeNewton(x1, xguess, x2, xacc, var, maxNRStep);
472 return optx;
473 }
474
475 /*****************************************************
476 Multi dimensional optimization with Newton Raphson
477 only applicable if 1st and 2nd derivatives are easy to compute
478 *****************************************************/
479
480 int matinv (double x[], int n, int m, double space[]);
481
minimizeNewtonMulti(double * x1,double * xguess,double * x2,double xacc,int N,int maxNRStep)482 double Optimization::minimizeNewtonMulti(double *x1, double *xguess, double *x2, double xacc, int N, int maxNRStep)
483 {
484 int i, step;
485 int N_2 = N*N;
486 double df[N_2], dx[N], dxold[N], f[N+1], space[N];
487 double temp, xh[N], xl[N], rts[N], rts_old[N];
488
489 for (i = 0; i < N; i++) {
490 rts[i] = xguess[i];
491 if (rts[i] < x1[i]) rts[i] = x1[i];
492 if (rts[i] > x2[i]) rts[i] = x2[i];
493 }
494 computeFuncDervMulti(rts, f, df);
495 double score = f[N];
496 // double score = targetFunk(rts-1);
497 // cout << "Newton step 0: " << score << endl;
498
499 bool stop = true;
500
501 for (i = 0; i < N; i++) {
502 if (!(fabs(f[i]) < xacc))
503 stop = false;
504 xl[i] = x1[i];
505 xh[i] = x2[i];
506 // if (f[i] < 0.0) {
507 // xl[i] = rts[i];
508 // xh[i] = x2[i];
509 // } else {
510 // xh[i] = rts[i];
511 // xl[i] = x1[i];
512 // }
513 dx[i] = dxold[i] = fabs(xh[i]-xl[i]);
514 }
515
516 if (stop) {
517 memcpy(xguess, rts, sizeof(double)*N);
518 return score;
519 }
520
521 for (step = 1; step <= maxNRStep; step++) {
522 memcpy(rts_old, rts, sizeof(double)*N);
523
524 // first compute inverse of hessian matrix
525 matinv(df, N, N, space);
526 // now take the product
527 for (i = 0; i < N; i++) {
528 space[i] = 0.0;
529 for (int k = 0; k < N; k++) {
530 space[i] += df[i*N+k] * f[k];
531 }
532 }
533
534 // identify alpha to be not out of bound
535 double alpha;
536 for (alpha = 1.0; alpha > xacc; alpha *= 0.5) {
537 bool ok_bound = true;
538 for (i = 0; i < N; i++)
539 if (rts[i]-alpha*space[i] < xl[i] || rts[i]-alpha*space[i] > xh[i]) {
540 ok_bound = false;
541 break;
542 }
543 if (ok_bound) break;
544 }
545
546 do {
547 stop = true;
548 for (i = 0; i < N; i++) {
549 dxold[i] = dx[i];
550 dx[i] = alpha*space[i];
551 temp=rts[i];
552 rts[i] -= dx[i];
553 if (fabs(dx[i]) >= xacc && (step != maxNRStep)) {
554 stop = false;
555 } else {
556 rts[i] = rts_old[i];
557 }
558 }
559 double new_score = targetFunk(rts-1);
560 if (new_score <= score) break;
561 // score increased, reduce step size by half
562 memcpy(rts, rts_old, sizeof(double)*N);
563 alpha *= 0.5;
564 } while (!stop);
565
566 if (stop) {
567 break;
568 }
569 computeFuncDervMulti(rts, f, df);
570 score = f[N];
571
572 // double score = targetFunk(rts-1);
573 // cout << "Newton step " << step << ": " << score << endl;
574 // for (i = 0; i < N; i++) cout << " " << rts[i];
575 // cout << endl;
576
577 stop = true;
578 for (i = 0; i < N; i++) {
579 if (!(fabs(f[i]) < xacc)) {
580 stop = false;
581 }
582 // if (f[i] < 0.0)
583 // xl[i] = rts[i];
584 // else
585 // xh[i] = rts[i];
586 }
587 if (stop)
588 break;
589 }
590
591 // copy returned x-value
592 memcpy(xguess, rts, sizeof(double)*N);
593
594 if (stop)
595 return score;
596
597 cout << "Maximum number of iterations exceeded in minimizeNewtonMulti" << endl;
598 return score;
599 }
600
601
602 /*****************************************************
603 Multi dimensional optimization with BFGS method
604 *****************************************************/
605
606 #define ALF 1.0e-4
607 #define TOLX 1.0e-7
608 //static double maxarg1,maxarg2;
609 //#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
610 // (maxarg1) : (maxarg2))
611
lnsrch(int n,double xold[],double fold,double g[],double p[],double x[],double * f,double stpmax,int * check,double lower[],double upper[])612 void Optimization::lnsrch(int n, double xold[], double fold, double g[], double p[], double x[],
613 double *f, double stpmax, int *check, double lower[], double upper[]) {
614 int i;
615 double a,alam,alam2=0,alamin,b,disc,f2=0,fold2=0,rhs1,rhs2,slope,sum,temp,
616 test,tmplam;
617
618 *check=0;
619 for (sum=0.0,i=1;i<=n;i++) sum += p[i]*p[i];
620 sum=sqrt(sum);
621 if (sum > stpmax)
622 for (i=1;i<=n;i++) p[i] *= stpmax/sum;
623 for (slope=0.0,i=1;i<=n;i++)
624 slope += g[i]*p[i];
625 test=0.0;
626 for (i=1;i<=n;i++) {
627 temp=fabs(p[i])/max(fabs(xold[i]),1.0);
628 if (temp > test) test=temp;
629 }
630 alamin=TOLX/test;
631 alam=1.0;
632 bool first_time = true;
633 for (;;) {
634 // for (; alam >= alamin;) {
635 // bool out_of_bound = false;
636 // for (i=1;i<=n;i++) {
637 // x[i]=xold[i]+alam*p[i];
638 // if (x[i] < lower[i] || x[i] > upper[i]) {
639 // out_of_bound = true;
640 // break;
641 // }
642 // }
643 // if (out_of_bound)
644 // alam *= 0.5;
645 // else
646 // break;
647 // }
648
649 for (i=1;i<=n;i++) x[i]=xold[i]+alam*p[i];
650
651 fixBound(x, lower, upper, n);
652 //checkRange(x);
653 *f=targetFunk(x);
654 if (alam < alamin) {
655 for (i=1;i<=n;i++) x[i]=xold[i];
656 *check=1;
657 return;
658 } else if (*f <= fold+ALF*alam*slope) return;
659 else {
660 if (first_time)
661 tmplam = -slope/(2.0*(*f-fold-slope));
662 else {
663 rhs1 = *f-fold-alam*slope;
664 rhs2=f2-fold2-alam2*slope;
665 a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
666 b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
667 if (a == 0.0) tmplam = -slope/(2.0*b);
668 else {
669 disc=b*b-3.0*a*slope;
670 if (disc<0.0) //nrerror("Roundoff problem in lnsrch.");
671 tmplam = 0.5 * alam;
672 else if (b <= 0.0) tmplam=(-b+sqrt(disc))/(3.0*a);
673 else tmplam = -slope/(b+sqrt(disc));
674 }
675 if (tmplam>0.5*alam)
676 tmplam=0.5*alam;
677 }
678 }
679 alam2=alam;
680 f2 = *f;
681 fold2=fold;
682 alam=max(tmplam,0.1*alam);
683 first_time = false;
684 }
685 }
686 #undef ALF
687 #undef TOLX
688
689
690 const int MAX_ITER = 3;
691 extern double random_double(int *rstream);
692
restartParameters(double guess[],int ndim,double lower[],double upper[],bool bound_check[],int iteration)693 bool Optimization::restartParameters(double guess[], int ndim, double lower[], double upper[], bool bound_check[], int iteration) {
694 int i;
695 bool restart = false;
696 if (iteration <= MAX_ITER) {
697 for (i = 1; i <= ndim; i++) {
698 if (bound_check[i]) {
699 if (fabs(guess[i]-lower[i]) < 1e-4 || fabs(guess[i]-upper[i]) < 1e-4) {
700 restart = true;
701 break;
702 } // if (fabs...
703 } // if bound_check
704 } // for
705 } // if iteration
706 if (restart) {
707 // TODO: Dominik observed that this warning is printed even during model
708 // testing. A quieter solution might be preferred.
709 cout << "Restart estimation at the boundary... " << std::endl;
710 for (i = 1; i <= ndim; i++) {
711 guess[i] = random_double() * (upper[i] - lower[i])/3 + lower[i];
712 }
713 }
714 return (restart);
715 }
716
minimizeMultiDimen(double guess[],int ndim,double lower[],double upper[],bool bound_check[],double gtol,double * hessian)717 double Optimization::minimizeMultiDimen(double guess[], int ndim, double lower[], double upper[], bool bound_check[], double gtol, double *hessian) {
718 int i, iter;
719 double fret, minf = 1e+12;
720 double *minx = new double [ndim+1];
721 int count = 0;
722 bool restart;
723 do {
724 dfpmin(guess, ndim, lower, upper, gtol, &iter, &fret, hessian);
725 if (fret < minf) {
726 minf = fret;
727 for (i = 1; i <= ndim; i++)
728 minx[i] = guess[i];
729 }
730 count++;
731 // restart the search if at the boundary
732 // it's likely to end at a local optimum at the boundary
733 restart = restartParameters(guess,ndim,lower,upper,bound_check,count);
734 } while (restart);
735 if (count > 1) {
736 for (i = 1; i <= ndim; i++)
737 guess[i] = minx[i];
738 fret = minf;
739 }
740 delete [] minx;
741
742 return fret;
743 }
744
745
746 #define ITMAX 200
747 #define SQR(a) ((a)*(a))
748 #define EPS 3.0e-8
749 #define TOLX (4*EPS)
750 #define STPMX 100.0
751
752 #define FREEALL free_vector(xi,1,n);free_vector(pnew,1,n); \
753 free_matrix(hessin,1,n,1,n);free_vector(hdg,1,n);free_vector(g,1,n); \
754 free_vector(dg,1,n);
755
756
757
dfpmin(double p[],int n,double lower[],double upper[],double gtol,int * iter,double * fret,double * hessian)758 void Optimization::dfpmin(double p[], int n, double lower[], double upper[], double gtol, int *iter, double *fret, double *hessian) {
759 int check,i,its,j;
760 double den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test;
761 double *dg,*g,*hdg,**hessin,*pnew,*xi;
762
763 dg=new_vector(1,n);
764 g=new_vector(1,n);
765 hdg=new_vector(1,n);
766 hessin=new_matrix(1,n,1,n);
767 pnew=new_vector(1,n);
768 xi=new_vector(1,n);
769 fp = derivativeFunk(p,g);
770 for (i=1;i<=n;i++) {
771 for (j=1;j<=n;j++) hessin[i][j]=0.0;
772 hessin[i][i]=1.0;
773 xi[i] = -g[i];
774 sum += p[i]*p[i];
775 }
776 if (hessian) {
777 for (i=1; i<=n; i++)
778 for (j=1; j<=n; j++)
779 hessin[i][j] = hessian[(i-1)*n+j-1];
780 }
781 //checkBound(p, xi, lower, upper, n);
782 //checkDirection(p, xi);
783
784 stpmax=STPMX*max(sqrt(sum),(double)n);
785 for (its=1;its<=ITMAX;its++) {
786 *iter=its;
787 lnsrch(n,p,fp,g,xi,pnew,fret,stpmax,&check, lower, upper);
788 fp = *fret;
789 for (i=1;i<=n;i++) {
790 xi[i]=pnew[i]-p[i];
791 p[i]=pnew[i];
792 }
793 test=0.0;
794 for (i=1;i<=n;i++) {
795 temp=fabs(xi[i])/max(fabs(p[i]),1.0);
796 if (temp > test) test=temp;
797 }
798 if (test < TOLX) {
799 FREEALL
800 return;
801 }
802 for (i=1;i<=n;i++) dg[i]=g[i];
803 derivativeFunk(p,g);
804 test=0.0;
805 den=max(fabs(*fret),1.0); // fix bug found by Tung, as also suggested by NR author
806 for (i=1;i<=n;i++) {
807 temp=fabs(g[i])*max(fabs(p[i]),1.0)/den;
808 if (temp > test) test=temp;
809 }
810 if (test < gtol) {
811 FREEALL
812 return;
813 }
814 for (i=1;i<=n;i++) dg[i]=g[i]-dg[i];
815 for (i=1;i<=n;i++) {
816 hdg[i]=0.0;
817 for (j=1;j<=n;j++) hdg[i] += hessin[i][j]*dg[j];
818 }
819 fac=fae=sumdg=sumxi=0.0;
820 for (i=1;i<=n;i++) {
821 fac += dg[i]*xi[i];
822 fae += dg[i]*hdg[i];
823 sumdg += SQR(dg[i]);
824 sumxi += SQR(xi[i]);
825 }
826 if (fac*fac > EPS*sumdg*sumxi)
827 {
828 fac=1.0/fac;
829 fad=1.0/fae;
830 for (i=1;i<=n;i++) dg[i]=fac*xi[i]-fad*hdg[i];
831 for (i=1;i<=n;i++) {
832 for (j=1;j<=n;j++) {
833 hessin[i][j] += fac*xi[i]*xi[j]
834 -fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
835 }
836 }
837 }
838 for (i=1;i<=n;i++) {
839 xi[i]=0.0;
840 for (j=1;j<=n;j++) xi[i] -= hessin[i][j]*g[j];
841 }
842 //checkBound(p, xi, lower, upper, n);
843 //checkDirection(p, xi);
844 //if (*iter > 200) cout << "iteration=" << *iter << endl;
845 }
846 // BQM: disable this message!
847 //nrerror("too many iterations in dfpmin");
848 FREEALL
849 }
850 #undef ITMAX
851 #undef SQR
852 #undef EPS
853 #undef TOLX
854 #undef STPMX
855 #undef FREEALL
856 #undef FMAX
857
858
859 /**
860 the approximated derivative function
861 @param x the input vector x
862 @param dfx the derivative at x
863 @return the function value at x
864 */
derivativeFunk(double x[],double dfx[])865 double Optimization::derivativeFunk(double x[], double dfx[]) {
866 /*
867 if (!checkRange(x))
868 return INFINITIVE;
869 */
870 int ndim = getNDim();
871 double *h = new double[ndim+1];
872 double temp;
873 int dim;
874 double fx = targetFunk(x);
875 for (dim = 1; dim <= ndim; dim++ ){
876 temp = x[dim];
877 h[dim] = ERROR_X * fabs(temp);
878 if (h[dim] == 0.0) h[dim] = ERROR_X;
879 x[dim] = temp + h[dim];
880 h[dim] = x[dim] - temp;
881 dfx[dim] = (targetFunk(x));
882 x[dim] = temp;
883 }
884 for (dim = 1; dim <= ndim; dim++ )
885 dfx[dim] = (dfx[dim] - fx) / h[dim];
886 delete [] h;
887 return fx;
888 }
889
890
891 /*#define NRANSI
892 #define ITMAX 100
893 #define CGOLD 0.3819660
894 #define ZEPS 1.0e-10
895 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
896 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
897
898 double Optimization::brent(double ax, double bx, double cx, double tol,
899 double *xmin)
900 {
901 int iter;
902 double a,b,d=0.0,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
903 double e=0.0;
904
905 a=(ax < cx ? ax : cx);
906 b=(ax > cx ? ax : cx);
907 x=w=v=bx;
908 fw=fv=fx=computeFunction(x);
909 for (iter=1;iter<=ITMAX;iter++) {
910 xm=0.5*(a+b);
911 tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
912 if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
913 *xmin=x;
914 return fx;
915 }
916 if (fabs(e) > tol1) {
917 r=(x-w)*(fx-fv);
918 q=(x-v)*(fx-fw);
919 p=(x-v)*q-(x-w)*r;
920 q=2.0*(q-r);
921 if (q > 0.0) p = -p;
922 q=fabs(q);
923 etemp=e;
924 e=d;
925 if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
926 d=CGOLD*(e=(x >= xm ? a-x : b-x));
927 else {
928 d=p/q;
929 u=x+d;
930 if (u-a < tol2 || b-u < tol2)
931 d=SIGN(tol1,xm-x);
932 }
933 } else {
934 d=CGOLD*(e=(x >= xm ? a-x : b-x));
935 }
936 u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
937 fu=computeFunction(u);
938 if (fu <= fx) {
939 if (u >= x) a=x; else b=x;
940 SHFT(v,w,x,u)
941 SHFT(fv,fw,fx,fu)
942 } else {
943 if (u < x) a=u; else b=u;
944 if (fu <= fw || w == x) {
945 v=w;
946 w=u;
947 fv=fw;
948 fw=fu;
949 } else if (fu <= fv || v == x || v == w) {
950 v=u;
951 fv=fu;
952 }
953 }
954 }
955 nrerror("Too many iterations in brent");
956 *xmin=x;
957 return fx;
958 }
959
960 #undef SIGN
961 #undef ITMAX
962 #undef CGOLD
963 #undef ZEPS
964 #undef SHFT
965 #undef NRANSI*/
966
967 /*#define JMAX 20
968
969 double Optimization::minimizeNewton(double xmin, double xguess, double xmax, double tolerance, double &f)
970 {
971 return rtsafe(xmin, xguess, xmax, tolerance, f);
972 //double fe;
973 //return minimizeOneDimen(xmin, rtn, xmax, tolerance, &f, &fe);
974
975 int j;
976 double df,ddf,dx,rtn,rtnold, fstart=0, fnew;
977
978 rtn=xguess;
979 if (rtn < xmin) rtn = xmin;
980 if (rtn > xmax) rtn = xmax;
981
982
983 for (j=1;j<=JMAX;j++) {
984 f = computeFuncDerv(rtn,df,ddf);
985 if (!isfinite(f))
986 return 0;
987 if (j == 1) fstart = f;
988 if (ddf == 0.0) break;
989 dx=(df/fabs(ddf));
990 if (fabs(dx) <= tolerance) break;
991
992 rtnold = rtn; rtn = rtn-dx;
993 if (rtn < xmin) rtn = xmin;
994 if (rtn > xmax) rtn = xmax;
995 dx = rtnold-rtn;
996
997 while (fabs(dx) > tolerance && (fnew = computeFunction(rtn)) > f + tolerance) {
998 dx /= 2;
999 rtn = rtnold - dx;
1000 }
1001 if (fabs(dx) <= tolerance) { rtn = rtnold; break; }
1002 }
1003 //if (j > JMAX)
1004 //nrerror("Maximum number of iterations exceeded in Newton-Raphson");
1005 if (f <= fstart && j <= JMAX && (j > 1 || xguess > xmin+tolerance))
1006 return rtn;
1007 // Newton does not work, turn to other method
1008 double fe;
1009 return minimizeOneDimen(xmin, xguess, xmax, tolerance, &f, &fe);
1010 }*/
1011
1012 /*
1013 double Optimization::minimizeNewton(double xmin, double xguess, double xmax, double tolerance, double &f)
1014 {
1015 int j;
1016 double df,ddf,dx,dxold,rtn,temp, fold, fstart;
1017 double xmin_orig = xmin, xmax_orig = xmax;
1018
1019 rtn=xguess;
1020 dx=dxold=(xmax-xmin);
1021 fstart = fold = f = computeFuncDerv(rtn,df,ddf);
1022
1023 for (j=1;j<=JMAX;j++) {
1024 if (ddf <= 0.0) break;
1025 if ((((rtn-xmax)*ddf-df) * ((rtn-xmin)*ddf-df) > 0) || // run out of range
1026 (fabs(2.0*df) > fabs(dxold*ddf)) // dx not decreasing fast enough
1027 ) // f even increase
1028 {
1029 dxold = dx;
1030 dx = 0.5*(xmax-xmin);
1031 rtn = xmin+dx;
1032 if (xmin == rtn) break;
1033 } else {
1034 dxold=dx;
1035 if (ddf == 0.0)
1036 nrerror("2nd derivative is zero");
1037 dx=df/ddf;
1038 temp=rtn;
1039 //if (f > fold) dx /= 2.0;
1040 //if (ddf < 0) dx = -dx;
1041 rtn -= dx;
1042 //if (rtn < xmin) rtn = xmin;
1043 //if (rtn > xmax) rtn = xmax;
1044 dx = temp - rtn;
1045 if (temp == rtn) break;
1046 }
1047 if (fabs(dx) < tolerance) break;
1048 fold = f;
1049 f = computeFuncDerv(rtn,df,ddf);
1050 //if (f > fold) break; // Does not decrease function, escape
1051 if (df < 0.0)
1052 xmin = rtn;
1053 else
1054 xmax = rtn;
1055 }
1056 if (j > JMAX)
1057 nrerror("Maximum number of iterations exceeded in Newton-Raphson");
1058 if (f <= fstart) return rtn;
1059 // Newton does not work (find a max instead of min), turn to other method
1060 double fe;
1061 return minimizeOneDimen(xmin_orig, xguess, xmax_orig, tolerance, &f, &fe);
1062 //return 0.0;
1063 }
1064 #undef JMAX*/
1065
1066
L_BFGS_B(int n,double * x,double * l,double * u,double pgtol,int maxit)1067 double Optimization::L_BFGS_B(int n, double* x, double* l, double* u, double pgtol, int maxit) {
1068 int i;
1069 double Fmin;
1070 int fail;
1071 int fncount;
1072 int grcount;
1073 char msg[100];
1074
1075 int m = 10; // number of BFGS updates retained in the "L-BFGS-B" method. It defaults to 5.
1076
1077 int *nbd; // 0: unbounded; 1: lower bounded; 2: both lower & upper; 3: upper bounded
1078 nbd = new int[n];
1079 for (i=0; i<n; i++)
1080 nbd[i] = 2;
1081
1082 double factr = 1e+7; // control the convergence of the "L-BFGS-B" method.
1083 // Convergence occurs when the reduction in the object is within this factor
1084 // of the machine tolerance.
1085 // Default is 1e7, that is a tolerance of about 1e-8
1086
1087 // double pgtol = 0; // helps control the convergence of the "L-BFGS-B" method.
1088 // pgtol = 0.0;
1089 // It is a tolerance on the projected gradient in the current search direction.
1090 // Default is zero, when the check is suppressed
1091
1092 int trace = 0; // non-negative integer.
1093 if (verbose_mode >= VB_MAX)
1094 trace = 1;
1095 // If positive, tracing information on the progress of the optimization is produced.
1096 // Higher values may produce more tracing information.
1097
1098 int nREPORT = 10; // The frequency of reports for the "L-BFGS-B" methods if "trace" is positive.
1099 // Defaults to every 10 iterations.
1100
1101 /*#ifdef USE_OLD_PARAM
1102 lbfgsb(n, m, x, l, u, nbd, &Fmin, fn, gr1, &fail, ex,
1103 factr, pgtol, &fncount, &grcount, maxit, msg, trace, nREPORT);
1104 #else*/
1105
1106 lbfgsb(n, m, x, l, u, nbd, &Fmin, &fail,
1107 factr, pgtol, &fncount, &grcount, maxit, msg, trace, nREPORT);
1108 //#endif
1109
1110 if (fail == 51 || fail == 52) {
1111 cout << msg << endl;
1112 }
1113
1114 delete[] nbd;
1115
1116 return Fmin;
1117 }
1118
lbfgsb(int n,int m,double * x,double * l,double * u,int * nbd,double * Fmin,int * fail,double factr,double pgtol,int * fncount,int * grcount,int maxit,char * msg,int trace,int nREPORT)1119 void Optimization::lbfgsb(int n, int m, double *x, double *l, double *u, int *nbd,
1120 double *Fmin, int *fail,
1121 double factr, double pgtol,
1122 int *fncount, int *grcount, int maxit, char *msg,
1123 int trace, int nREPORT)
1124 {
1125 char task[60];
1126 double f, *g, dsave[29], *wa;
1127 int tr = -1, iter = 0, *iwa, isave[44], lsave[4];
1128
1129 /* shut up gcc -Wall in 4.6.x */
1130
1131 for(int i = 0; i < 4; i++) lsave[i] = 0;
1132
1133 if(n == 0) { /* not handled in setulb */
1134 *fncount = 1;
1135 *grcount = 0;
1136 *Fmin = optimFunc(n, u);
1137 strcpy(msg, "NOTHING TO DO");
1138 *fail = 0;
1139 return;
1140 }
1141 if (nREPORT <= 0) {
1142 cerr << "REPORT must be > 0 (method = \"L-BFGS-B\")" << endl;
1143 exit(1);
1144 }
1145 switch(trace) {
1146 case 2: tr = 0; break;
1147 case 3: tr = nREPORT; break;
1148 case 4: tr = 99; break;
1149 case 5: tr = 100; break;
1150 case 6: tr = 101; break;
1151 default: tr = -1; break;
1152 }
1153
1154 *fail = 0;
1155 g = (double*) malloc (n * sizeof(double));
1156 /* this needs to be zeroed for snd in mainlb to be zeroed */
1157 wa = (double *) malloc((2*m*n+4*n+11*m*m+8*m) * sizeof(double));
1158 iwa = (int *) malloc(3*n * sizeof(int));
1159 strcpy(task, "START");
1160 while(1) {
1161 /* Main workhorse setulb() from ../appl/lbfgsb.c : */
1162 setulb(n, m, x, l, u, nbd, &f, g, factr, &pgtol, wa, iwa, task,
1163 tr, lsave, isave, dsave);
1164 /* Rprintf("in lbfgsb - %s\n", task);*/
1165 if (strncmp(task, "FG", 2) == 0) {
1166 f = optimGradient(n, x, g);
1167 if (!isfinite(f)) {
1168 cerr << "L-BFGS-B needs finite values of 'fn'" << endl;
1169 exit(1);
1170 }
1171
1172 } else if (strncmp(task, "NEW_X", 5) == 0) {
1173 iter++;
1174 if(trace == 1 && (iter % nREPORT == 0)) {
1175 cout << "iter " << iter << " value " << f << endl;
1176 }
1177 if (iter > maxit) {
1178 *fail = 1;
1179 break;
1180 }
1181 } else if (strncmp(task, "WARN", 4) == 0) {
1182 *fail = 51;
1183 break;
1184 } else if (strncmp(task, "CONV", 4) == 0) {
1185 break;
1186 } else if (strncmp(task, "ERROR", 5) == 0) {
1187 *fail = 52;
1188 break;
1189 } else { /* some other condition that is not supposed to happen */
1190 *fail = 52;
1191 break;
1192 }
1193 }
1194 *Fmin = f;
1195 *fncount = *grcount = isave[33];
1196 if (trace) {
1197 cout << "final value " << *Fmin << endl;
1198 if (iter < maxit && *fail == 0)
1199 cout << "converged" << endl;
1200 else
1201 cout << "stopped after " << iter << " iterations\n";
1202 }
1203 strcpy(msg, task);
1204 free(g);
1205 free(wa);
1206 free(iwa);
1207 }
1208
optimFunc(int nvar,double * vars)1209 double Optimization::optimFunc(int nvar, double *vars) {
1210 return targetFunk(vars-1);
1211 }
1212
1213
optimGradient(int nvar,double * x,double * dfx)1214 double Optimization::optimGradient(int nvar, double *x, double *dfx) {
1215 return derivativeFunk(x-1, dfx-1);
1216 // const double ERRORX = 1e-5;
1217 // double fx = optimFunc(nvar, x);
1218 // double h, temp;
1219 // for (int dim = 0; dim <= nvar; dim++ ){
1220 // temp = x[dim];
1221 // h = ERRORX * fabs(temp);
1222 // if (h == 0.0) h = ERRORX;
1223 // x[dim] = temp + h;
1224 // h = x[dim] - temp;
1225 // dfx[dim] = (optimFunc(nvar, x) - fx) / h;
1226 // x[dim] = temp;
1227 // }
1228 // return fx;
1229 }
1230