--- imach/src/imach.c 2014/09/15 20:41:41 1.161 +++ imach/src/imach.c 2014/09/25 11:43:39 1.162 @@ -1,6 +1,14 @@ -/* $Id: imach.c,v 1.161 2014/09/15 20:41:41 brouard Exp $ +/* $Id: imach.c,v 1.162 2014/09/25 11:43:39 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.162 2014/09/25 11:43:39 brouard + Summary: temporary backup 0.99! + + 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 @@ -515,6 +523,13 @@ #include #endif +#ifdef NLOPT +#include +typedef struct { + double (* function)(double [] ); +} myfunc_data ; +#endif + /* #include */ /* #define _(String) gettext (String) */ @@ -550,11 +565,11 @@ #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.161 2014/09/15 20:41:41 brouard Exp $ */ +/* $Id: imach.c,v 1.162 2014/09/25 11:43:39 brouard Exp $ */ /* $State: Exp $ */ -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 fullversion[]="$Revision: 1.161 $ $Date: 2014/09/15 20:41:41 $"; +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: 1.162 $ $Date: 2014/09/25 11:43:39 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -584,6 +599,7 @@ int **mw; /* mw[mi][i] is number of the 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 */ @@ -1093,6 +1109,19 @@ char *subdirf3(char fileres[], char *pre 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; @@ -1131,7 +1160,7 @@ double brent(double ax, double bx, doubl /* 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)))) { */ @@ -1201,21 +1230,21 @@ void mnbrak(double *ax, double *bx, doub } *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 { @@ -1228,7 +1257,11 @@ void mnbrak(double *ax, double *bx, doub } /*************** 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 []); @@ -1254,8 +1287,8 @@ void linmin(double p[], double xi[], int } 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); @@ -1268,20 +1301,16 @@ void linmin(double p[], double xi[], int 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 [])) { @@ -1322,17 +1351,15 @@ void powell(double p[], double **xi, int 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') @@ -1364,13 +1391,13 @@ void powell(double p[], double **xi, int 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; @@ -1410,14 +1437,17 @@ void powell(double p[], double **xi, int } 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); @@ -1447,8 +1477,8 @@ void powell(double p[], double **xi, int printf("\n"); fprintf(ficlog,"\n"); #endif - } - } + } /* end of t negative */ + } /* end if (fptt < fp) */ } } @@ -1678,6 +1708,25 @@ double ***hpxij(double ***po, int nhstep 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) @@ -1696,6 +1745,9 @@ double func( double *x) /*for(i=1;ifunction = 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)); } @@ -6019,7 +6110,7 @@ Interval (in months) between two waves: #ifdef GSL printf("GSL optimization\n"); fprintf(ficlog,"Powell\n"); -#elsedef +#else printf("Powell\n"); fprintf(ficlog,"Powell\n"); #endif strcpy(filerespow,"pow-mort"); @@ -6030,7 +6121,7 @@ Interval (in months) between two waves: } #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++)