Summary: Comments on Powell added

Author:

diff --git a/src/imach.c b/src/imach.c
index c20f7e63664b9815b1991f9949605b32858e6ad4..235b7f1f210112661971b1147c3d07a2e547b8bf 100644 (file)
--- 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 <math.h>
 #include <stdio.h>
@@ -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);