/* $Id$
$State$
$Log$
+ Revision 1.1 2014/09/16 11:06:58 brouard
+ Summary: With some code (wrong) for nlopt
+
+ Author:
+
+ Revision 1.161 2014/09/15 20:41:41 brouard
+ Summary: Problem with macro SQR on Intel compiler
+
Revision 1.160 2014/09/02 09:24:05 brouard
*** empty log message ***
#include <gsl/gsl_multimin.h>
#endif
+#ifdef NLOPT
+#include <nlopt.h>
+typedef struct {
+ double (* function)(double [] );
+} myfunc_data ;
+#endif
+
/* #include <libintl.h> */
/* #define _(String) gettext (String) */
/* $Id$ */
/* $State$ */
-char version[]="Imach version 0.98nX, August 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)";
+char version[]="Imach version 0.99, September 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)";
char fullversion[]="$Revision$ $Date$";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */
int **bh; /* bh[mi][i] is the bias (+ or -) for this individual if the delay between
* wave mi and wave mi+1 is not an exact multiple of stepm. */
+int countcallfunc=0; /* Count the number of calls to func */
double jmean=1; /* Mean space between 2 waves */
double **matprod2(); /* test */
double **oldm, **newm, **savm; /* Working pointers to matrices */
return tmpout;
}
+char *asc_diff_time(long time_sec, char ascdiff[])
+{
+ long sec_left, days, hours, minutes;
+ days = (time_sec) / (60*60*24);
+ sec_left = (time_sec) % (60*60*24);
+ hours = (sec_left) / (60*60) ;
+ sec_left = (sec_left) %(60*60);
+ minutes = (sec_left) /60;
+ sec_left = (sec_left) % (60);
+ sprintf(ascdiff,"%ld day(s) %ld hour(s) %ld minute(s) %ld second(s)",days, hours, minutes, sec_left);
+ return ascdiff;
+}
+
/***************** f1dim *************************/
extern int ncom;
extern double *pcom,*xicom;
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret)))*/
printf(".");fflush(stdout);
fprintf(ficlog,".");fflush(ficlog);
-#ifdef DEBUG
+#ifdef DEBUGBRENT
printf("br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);
fprintf(ficlog,"br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);
/* if ((fabs(x-xm) <= (tol2-0.5*(b-a)))||(2.0*fabs(fu-ftemp) <= ftol*1.e-2*(fabs(fu)+fabs(ftemp)))) { */
}
*cx=(*bx)+GOLD*(*bx-*ax);
*fc=(*func)(*cx);
- while (*fb > *fc) {
+ while (*fb > *fc) { /* Declining fa, fb, fc */
r=(*bx-*ax)*(*fb-*fc);
q=(*bx-*cx)*(*fb-*fa);
u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
- (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));
- ulim=(*bx)+GLIMIT*(*cx-*bx);
- if ((*bx-u)*(u-*cx) > 0.0) {
+ (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); /* Minimum abscisse of a parabolic estimated from (a,fa), (b,fb) and (c,fc). */
+ ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscisse where function can be evaluated */
+ if ((*bx-u)*(u-*cx) > 0.0) { /* if u between b and c */
fu=(*func)(u);
- } else if ((*cx-u)*(u-ulim) > 0.0) {
+ } else if ((*cx-u)*(u-ulim) > 0.0) { /* u is after c but before ulim */
fu=(*func)(u);
if (fu < *fc) {
SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
SHFT(*fb,*fc,fu,(*func)(u))
}
- } else if ((u-ulim)*(ulim-*cx) >= 0.0) {
+ } else if ((u-ulim)*(ulim-*cx) >= 0.0) { /* u outside ulim (verifying that ulim is beyond c) */
u=ulim;
fu=(*func)(u);
} else {
}
/*************** linmin ************************/
-
+/* Given an n -dimensional point p[1..n] and an n -dimensional direction xi[1..n] , moves and
+resets p to where the function func(p) takes on a minimum along the direction xi from p ,
+and replaces xi by the actual vector displacement that p was moved. Also returns as fret
+the value of func at the returned location p . This is actually all accomplished by calling the
+routines mnbrak and brent .*/
int ncom;
double *pcom,*xicom;
double (*nrfunc)(double []);
}
ax=0.0;
xx=1.0;
- mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);
- *fret=brent(ax,xx,bx,f1dim,TOL,&xmin);
+ mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Find a bracket a,x,b in direction n=xi ie xicom */
+ *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Find a minimum P+lambda n in that direction (lambdamin), with TOL between abscisses */
#ifdef DEBUG
printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
free_vector(pcom,1,n);
}
-char *asc_diff_time(long time_sec, char ascdiff[])
-{
- long sec_left, days, hours, minutes;
- days = (time_sec) / (60*60*24);
- sec_left = (time_sec) % (60*60*24);
- hours = (sec_left) / (60*60) ;
- sec_left = (sec_left) %(60*60);
- minutes = (sec_left) /60;
- sec_left = (sec_left) % (60);
- sprintf(ascdiff,"%ld day(s) %ld hour(s) %ld minute(s) %ld second(s)",days, hours, minutes, sec_left);
- return ascdiff;
-}
/*************** powell ************************/
+/*
+Minimization of a function func of n variables. Input consists of an initial starting point
+p[1..n] ; an initial matrix xi[1..n][1..n] , whose columns contain the initial set of di-
+rections (usually the n unit vectors); and ftol , the fractional tolerance in the function value
+such that failure to decrease by more than this amount on one iteration signals doneness. On
+output, p is set to the best point found, xi is the then-current direction set, fret is the returned
+function value at p , and iter is the number of iterations taken. The routine linmin is used.
+ */
void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret,
double (*func)(double []))
{
if(*iter <=3){
tml = *localtime(&rcurr_time);
strcpy(strcurr,asctime(&tml));
-/* asctime_r(&tm,strcurr); */
rforecast_time=rcurr_time;
itmp = strlen(strcurr);
if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */
strcurr[itmp-1]='\0';
- printf("\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
+ printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
for(niterf=10;niterf<=30;niterf+=10){
rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);
forecast_time = *localtime(&rforecast_time);
-/* asctime_r(&tmf,strfor); */
strcpy(strfor,asctime(&forecast_time));
itmp = strlen(strfor);
if(strfor[itmp-1]=='\n')
fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
}
for(j=1;j<=n;j++) {
- printf(" p=%.12e",p[j]);
- fprintf(ficlog," p=%.12e",p[j]);
+ printf(" p(%d)=%.12e",j,p[j]);
+ fprintf(ficlog," p(%d)=%.12e",j,p[j]);
}
printf("\n");
fprintf(ficlog,"\n");
#endif
- }
+ } /* end i */
if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
#ifdef DEBUG
int k[2],l;
}
fptt=(*func)(ptt);
if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */
- /* x1 f1=fp x2 f2=*fret x3 f3=fptt, xm fm */
+ /* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */
/* From x1 (P0) distance of x2 is at h and x3 is 2h */
- /* Let f"(x2) be the 2nd derivative equal everywhere. Then the parabolic through (x1,f1), (x2,f2) and (x3,f3)
- will reach at f3 = fm + h^2/2 f''m ; f" = (f1 -2f2 +f3 ) / h**2 */
+ /* Let f"(x2) be the 2nd derivative equal everywhere. */
+ /* Then the parabolic through (x1,f1), (x2,f2) and (x3,f3) */
+ /* will reach at f3 = fm + h^2/2 f"m ; f" = (f1 -2f2 +f3 ) / h**2 */
/* f1-f3 = delta(2h) = 2 h**2 f'' = 2(f1- 2f2 +f3) */
/* Thus we compare delta(2h) with observed f1-f3 */
- /* or best gain on one ancient line 'del' with total gain f1-f2 = f1 - f2 - 'del' with del */
+ /* or best gain on one ancient line 'del' with total */
+ /* gain f1-f2 = f1 - f2 - 'del' with del */
/* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */
+
t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del);
t= t- del*SQR(fp-fptt);
printf("t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t);
printf("\n");
fprintf(ficlog,"\n");
#endif
- }
- }
+ } /* end of t negative */
+ } /* end if (fptt < fp) */
}
}
return po;
}
+#ifdef NLOPT
+ double myfunc(unsigned n, const double *p1, double *grad, void *pd){
+ double fret;
+ double *xt;
+ int j;
+ myfunc_data *d2 = (myfunc_data *) pd;
+/* xt = (p1-1); */
+ xt=vector(1,n);
+ for (j=1;j<=n;j++) xt[j]=p1[j-1]; /* xt[1]=p1[0] */
+
+ fret=(d2->function)(xt); /* p xt[1]@8 is fine */
+ /* fret=(*func)(xt); /\* p xt[1]@8 is fine *\/ */
+ printf("Function = %.12lf ",fret);
+ for (j=1;j<=n;j++) printf(" %d %.8lf", j, xt[j]);
+ printf("\n");
+ free_vector(xt,1,n);
+ return fret;
+}
+#endif
/*************** log-likelihood *************/
double func( double *x)
/*for(i=1;i<imx;i++)
printf(" %d\n",s[4][i]);
*/
+
+ ++countcallfunc;
+
cov[1]=1.;
for(k=1; k<=nlstate; k++) ll[k]=0.;
double fret;
double fretone; /* Only one call to likelihood */
/* char filerespow[FILENAMELENGTH];*/
+
+#ifdef NLOPT
+ int creturn;
+ nlopt_opt opt;
+ /* double lb[9] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL }; /\* lower bounds *\/ */
+ double *lb;
+ double minf; /* the minimum objective value, upon return */
+ double * p1; /* Shifted parameters from 0 instead of 1 */
+ myfunc_data dinst, *d = &dinst;
+#endif
+
+
xi=matrix(1,npar,1,npar);
for (i=1;i<=npar;i++)
for (j=1;j<=npar;j++)
for(j=1;j<=nlstate+ndeath;j++)
if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
fprintf(ficrespow,"\n");
-
+#ifdef POWELL
powell(p,xi,npar,ftol,&iter,&fret,func);
+#endif
+#ifdef NLOPT
+#ifdef NEWUOA
+ opt = nlopt_create(NLOPT_LN_NEWUOA,npar);
+#else
+ opt = nlopt_create(NLOPT_LN_BOBYQA,npar);
+#endif
+ lb=vector(0,npar-1);
+ for (i=0;i<npar;i++) lb[i]= -HUGE_VAL;
+ nlopt_set_lower_bounds(opt, lb);
+ nlopt_set_initial_step1(opt, 0.1);
+
+ p1= (p+1); /* p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] */
+ d->function = func;
+ printf(" Func %.12lf \n",myfunc(npar,p1,NULL,d));
+ nlopt_set_min_objective(opt, myfunc, d);
+ nlopt_set_xtol_rel(opt, ftol);
+ if ((creturn=nlopt_optimize(opt, p1, &minf)) < 0) {
+ printf("nlopt failed! %d\n",creturn);
+ }
+ else {
+ printf("found minimum after %d evaluations (NLOPT=%d)\n", countcallfunc ,NLOPT);
+ printf("found minimum at f(%g,%g) = %0.10g\n", p[0], p[1], minf);
+ iter=1; /* not equal */
+ }
+ nlopt_destroy(opt);
+#endif
free_matrix(xi,1,npar,1,npar);
fclose(ficrespow);
- printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
- fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
- fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
+ printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
+ fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
+ fprintf(ficres,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
}
#ifdef GSL
printf("GSL optimization\n"); fprintf(ficlog,"Powell\n");
-#elsedef
+#else
printf("Powell\n"); fprintf(ficlog,"Powell\n");
#endif
strcpy(filerespow,"pow-mort");
}
#ifdef GSL
fprintf(ficrespow,"# GSL optimization\n# iter -2*LL");
-#elsedef
+#else
fprintf(ficrespow,"# Powell\n# iter -2*LL");
#endif
/* for (i=1;i<=nlstate;i++)
git://henry.ined.fr / .git / commitdiff