--- imach/src/imach.c	2015/04/29 09:11:15	1.187
+++ imach/src/imach.c	2015/04/30 08:27:53	1.188
@@ -1,6 +1,9 @@
-/* $Id: imach.c,v 1.187 2015/04/29 09:11:15 brouard Exp $
+/* $Id: imach.c,v 1.188 2015/04/30 08:27:53 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.188  2015/04/30 08:27:53  brouard
+  *** empty log message ***
+
   Revision 1.187  2015/04/29 09:11:15  brouard
   *** empty log message ***
 
@@ -675,11 +678,11 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.187 2015/04/29 09:11:15 brouard Exp $ */
+/* $Id: imach.c,v 1.188 2015/04/30 08:27:53 brouard Exp $ */
 /* $State: Exp $ */
 
 char version[]="Imach version 0.98q1, April 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.187 $ $Date: 2015/04/29 09:11:15 $"; 
+char fullversion[]="$Revision: 1.188 $ $Date: 2015/04/30 08:27:53 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1565,13 +1568,14 @@ void linmin(double p[], double xi[], int
 #endif
   printf("linmin end ");
   for (j=1;j<=n;j++) { 
-    printf(" before xi[%d]=%12.8f", j,xi[j]);
+    /* printf(" before xi[%d]=%12.8f", j,xi[j]); */
     xi[j] *= xmin; /* xi rescaled by xmin: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */
     if(xxs <1.0)
       printf(" after xi[%d]=%12.8f, xmin=%12.8f, ax=%12.8f, xx=%12.8f, bx=%12.8f, xxs=%12.8f", j,xi[j], xmin, ax, xx, bx,xxs );
     p[j] += xi[j]; /* Parameters values are updated accordingly */
   } 
   printf("\n");
+  printf("Comparing last *frec(xmin)=%12.8f from Brent and frec(0.)=%12.8f \n", (*func)(p));
   free_vector(xicom,1,n); 
   free_vector(pcom,1,n); 
 } 
@@ -1653,13 +1657,13 @@ void powell(double p[], double **xi, int
 #endif
 	  printf("%d",i);fflush(stdout); /* print direction (parameter) i */
       fprintf(ficlog,"%d",i);fflush(ficlog);
-      linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input. Outputs are fret(new point p) p is updated and xit rescaled */
-      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.
-				      */
+      linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input.*/
+				    /* Outputs are fret(new point p) p is updated and xit rescaled */
+      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; 
       } 
@@ -1680,9 +1684,21 @@ void powell(double p[], double **xi, int
 #endif
     } /* end loop on each direction i */
     /* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */ 
-    /* But p and xit have been updated at the end of linmin and do not produce *fret any more! */
+    /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit  */
     /* New value of last point Pn is not computed, P(n-1) */
     if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */
+      /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */
+      /* By adding age*age in a model, the new -2LL should be lower and the difference follows a */
+      /* a chisquare statistics with 1 degree. To be significant at the 95% level, it should have */
+      /* decreased of more than 3.84  */
+      /* By adding age*age and V1*age the gain (-2LL) should be more than 5.99 (ddl=2) */
+      /* By using V1+V2+V3, the gain should be  7.82, compared with basic 1+age. */
+      /* By adding 10 parameters more the gain should be 18.31 */
+
+      /* Starting the program with initial values given by a former maximization will simply change */
+      /* the scales of the directions and the directions, because the are reset to canonical directions */
+      /* Thus the first calls to linmin will give new points and better maximizations until fp-(*fret) is */
+      /* under the tolerance value. If the tolerance is very small 1.e-9, it could last long.  */
 #ifdef DEBUG
       int k[2],l;
       k[0]=1;