From: N. Brouard Date: Wed, 11 Feb 2015 23:22:24 +0000 (+0000) Subject: Summary: Comments on Powell added X-Git-Tag: imach-099s7~378 X-Git-Url: https://henry.ined.fr/git/?a=commitdiff_plain;h=e826e3ef116052b1f8180c67a083d59b08f266df;p=.git Summary: Comments on Powell added Author: --- diff --git a/src/imach.c b/src/imach.c index c20f7e6..235b7f1 100644 --- a/src/imach.c +++ b/src/imach.c @@ -1,6 +1,9 @@ /* $Id$ $State$ $Log$ + Revision 1.180 2015/02/11 17:33:45 brouard + Summary: Finishing move from main to function (hpijx and prevalence_limit) + Revision 1.179 2015/01/04 09:57:06 brouard Summary: back to OS/X @@ -553,6 +556,7 @@ */ #define POWELL /* Instead of NLOPT */ +#define POWELLDIRECT /* Directest to decide new direction instead of Powell test */ #include #include @@ -1401,6 +1405,7 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret double (*func)(double [])); int i,ibig,j; double del,t,*pt,*ptt,*xit; + double directest; double fp,fptt; double *xits; int niterf, itmp; @@ -1461,7 +1466,12 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret 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; } @@ -1513,29 +1523,27 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret 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), @@ -1543,14 +1551,18 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret 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);