From 75ffc90d2dd7542df5c19d51beb3587ec635ae91 Mon Sep 17 00:00:00 2001
From: "N. Brouard" <brouard@ined.fr>
Date: Thu, 30 Apr 2015 08:27:53 +0000
Subject: [PATCH] *** empty log message ***

---
 src/imach.c | 34 +++++++++++++++++++++++++---------
 1 file changed, 25 insertions(+), 9 deletions(-)

diff --git a/src/imach.c b/src/imach.c
index bd232dd..2646c25 100644
--- a/src/imach.c
+++ b/src/imach.c
@@ -1,6 +1,9 @@
 /* $Id$
   $State$
   $Log$
+  Revision 1.187  2015/04/29 09:11:15  brouard
+  *** empty log message ***
+
   Revision 1.186  2015/04/23 12:01:52  brouard
   Summary: V1*age is working now, version 0.98q1
 
@@ -1562,13 +1565,14 @@ void linmin(double p[], double xi[], int n, double *fret,double (*func)(double [
 #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); 
 } 
@@ -1650,13 +1654,13 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
 #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; 
       } 
@@ -1677,9 +1681,21 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
 #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;
-- 
2.43.0