--- imach/src/imach.c 2015/02/11 17:33:45 1.180 +++ imach/src/imach.c 2015/02/11 23:22:24 1.181 @@ -1,6 +1,11 @@ -/* $Id: imach.c,v 1.180 2015/02/11 17:33:45 brouard Exp $ +/* $Id: imach.c,v 1.181 2015/02/11 23:22:24 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.181 2015/02/11 23:22:24 brouard + Summary: Comments on Powell added + + Author: + Revision 1.180 2015/02/11 17:33:45 brouard Summary: Finishing move from main to function (hpijx and prevalence_limit) @@ -556,6 +561,7 @@ */ #define POWELL /* Instead of NLOPT */ +#define POWELLDIRECT /* Directest to decide new direction instead of Powell test */ #include #include @@ -638,11 +644,11 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.180 2015/02/11 17:33:45 brouard Exp $ */ +/* $Id: imach.c,v 1.181 2015/02/11 23:22:24 brouard Exp $ */ /* $State: Exp $ */ char version[]="Imach version 0.98p, Février 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; -char fullversion[]="$Revision: 1.180 $ $Date: 2015/02/11 17:33:45 $"; +char fullversion[]="$Revision: 1.181 $ $Date: 2015/02/11 23:22:24 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1404,6 +1410,7 @@ void powell(double p[], double **xi, int double (*func)(double [])); int i,ibig,j; double del,t,*pt,*ptt,*xit; + double directest; double fp,fptt; double *xits; int niterf, itmp; @@ -1464,7 +1471,12 @@ void powell(double p[], double **xi, int printf("%d",i);fflush(stdout); fprintf(ficlog,"%d",i);fflush(ficlog); linmin(p,xit,n,fret,func); - if (fabs(fptt-(*fret)) > del) { + if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions + because that direction will be replaced unless the gain del is small + in comparison with the 'probable' gain, mu^2, with the last average direction. + Unless the n directions are conjugate some gain in the determinant may be obtained + with the new direction. + */ del=fabs(fptt-(*fret)); ibig=i; } @@ -1516,29 +1528,27 @@ void powell(double p[], double **xi, int return; } if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); - for (j=1;j<=n;j++) { /* Computes an extrapolated point */ + for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */ ptt[j]=2.0*p[j]-pt[j]; xit[j]=p[j]-pt[j]; pt[j]=p[j]; } - fptt=(*func)(ptt); + fptt=(*func)(ptt); /* f_3 */ 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) */ /* 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 */ - /* 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 */ + /* Conditional for using this new direction is that mu^2 = (f1-2f2+f3)^2 /2 < 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); - fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t); + directest = SQR(fp-2.0*(*fret)+fptt) - 2.0 * del; /* If del was big enough we change it for a new direction */ #ifdef DEBUG + printf("t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest); + fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest); printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt), (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt)); fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt), @@ -1546,14 +1556,18 @@ void powell(double p[], double **xi, int printf("tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t); fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t); #endif - if (t < 0.0) { /* Then we use it for last direction */ - linmin(p,xit,n,fret,func); /* computes mean on the extrapolated direction.*/ +#ifdef POWELLDIRECT + if (directest < 0.0) { /* Then we use it for new direction */ +#else + if (t < 0.0) { /* Then we use it for new direction */ +#endif + linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction.*/ for (j=1;j<=n;j++) { - xi[j][ibig]=xi[j][n]; /* Replace the direction with biggest decrease by n */ - xi[j][n]=xit[j]; /* and nth direction by the extrapolated */ + xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */ + xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ } - printf("Gaining to use average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); - fprintf(ficlog,"Gaining to use average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); + printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); + fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); #ifdef DEBUG printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);