*** empty log message ***

diff --git a/src/imachprax.c b/src/imachprax.c
index e5da29d460d3a7db2c8467550552f6f3b6f2fb33..b20b69f32fec186f30bd420ee403f988d02641fe 100644 (file)
--- a/src/imachprax.c
+++ b/src/imachprax.c
@@ -1,6 +1,9 @@
 /* $Id$
   $State$
   $Log$
+  Revision 1.2  2023/06/22 11:22:40  brouard
+  Summary: with svd but not working yet
+
   Revision 1.353  2023/05/08 18:48:22  brouard
   *** empty log message ***
 
@@ -1278,6 +1281,7 @@ Important routines
 /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
 /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
 /* #define FLATSUP  *//* Suppresses directions where likelihood is flat */
+/* #define POWELLORIGINCONJUGATE  /\* Don't use conjugate but biggest decrease if valuable *\/ */
 
 #include <math.h>
 #include <stdio.h>
@@ -2996,7 +3000,9 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
       printf("%d",i);fflush(stdout); /* print direction (parameter) i */
       fprintf(ficlog,"%d",i);fflush(ficlog);
 #ifdef LINMINORIGINAL
-      linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input.*/
+      linmin(p,xit,n,fret,func); /* New point i minimizing in direction i has coordinates p[j].*/
+      /* xit[j] gives the n coordinates of direction i as input.*/
+      /* *fret gives the maximum value on direction xit */
 #else
       linmin(p,xit,n,fret,func,&flat); /* Point p[n]. xit[n] has been loaded for direction i as input.*/
                        flatdir[i]=flat; /* Function is vanishing in that direction i */
@@ -3026,9 +3032,8 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
       fprintf(ficlog,"\n");
 #endif
     } /* end loop on each direction i */
-    /* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */ 
+    /* Convergence test will use last linmin estimation (fret) and compare to former iteration (fp) */ 
     /* 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) */
     for(j=1;j<=n;j++) {
       if(flatdir[j] >0){
         printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
@@ -3166,10 +3171,17 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
          }
        }
 #endif
+#ifdef POWELLORIGINCONJUGATE
        for (j=1;j<=n;j++) { 
          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 */
        }
+#else
+       for (j=1;j<=n-1;j++) { 
+         xi[j][1]=xi[j][j+1]; /* Standard method of conjugate directions */
+         xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
+       }
+#endif
 #ifdef LINMINORIGINAL
 #else
        for (j=1, flatd=0;j<=n;j++) {
@@ -3198,6 +3210,15 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
 #endif
        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);
+  /* The minimization in direction $\xi_1$ gives $P_1$. From $P_1$ minimization in direction $\xi_2$ gives */
+  /* $P_2$. Minimization of line $P_2-P_1$ gives new starting point $P^{(1)}_0$ and direction $\xi_1$ is dropped and replaced by second */
+  /* direction $\xi_1^{(1)}=\xi_2$. Also second direction is replaced by new direction $\xi^{(1)}_2=P_2-P_0$. */
+
+  /* At the second iteration, starting from $P_0^{(1)}$, minimization amongst $\xi^{(1)}_1$ gives point $P^{(1)}_1$. */
+  /* Minimization amongst $\xi^{(1)}_2=(P_2-P_0)$ gives point $P^{(1)}_2$.  As $P^{(2)}_1$ and */
+  /* $P^{(1)}_0$ are minimizing in the same direction $P^{(1)}_2 - P^{(1)}_1= P_2-P_0$, directions $P^{(1)}_2-P^{(1)}_0$ */
+  /* and $P_2-P_0$ (parallel to $\xi$ and $\xi^c$) are conjugate.  } */
+
        
 #ifdef DEBUG
        printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
@@ -4808,7 +4829,7 @@ double funcone( double *x)
        * 3 ncovta=15    +age*V3*V2+age*V2+agev3+ageV4 +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4
        * 3 TvarAVVA[1]@15= itva 3 2    2      3    4        6       7        6 3         7 3         6 4         7 4 
        * 3 ncovta             1 2      3      4    5        6       7        8 9       10 11       12 13        14 15
-       * 
+       *?TvarAVVAind[1]@15= V3 is in k=2 1 1  2    3        4       5        4,2         5,2,      4,3           5 3}TvarVVAind[]
        * TvarAVVAind[1]@15= V3 is in k=6 6 12  13   14      15      16       18 18       19,19,     20,20        21,21}TvarVVAind[]
        * 3 ncovvta=10     +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4
        * 3 we want to compute =cotvar[mw[mi][i]][TvarVVA[ncovva]][i] at position TvarVVAind[ncovva]
@@ -4822,7 +4843,8 @@ double funcone( double *x)
        *                   2, 3, 4, 6, 7,
        *                     6, 8, 9, 10, 11}
        * TvarFind[itv]                        0      0       0
-       * FixedV[itv]                          1      1       1  0      1 0       1 0       1 0      1 0     1 0
+       * FixedV[itv]                          1      1       1  0      1 0       1 0       1 0       0
+       *? FixedV[itv]                          1      1       1  0      1 0       1 0       1 0      1 0     1 0
        * Tvar[TvarFind[ncovf]]=[1]=2 [2]=3 [4]=4
        * Tvar[TvarFind[itv]]                [0]=?      ?ncovv 1 à ncovvt]
        *   Not a fixed cotvar[mw][itv][i]     6       7      6  2      7, 2,     6, 3,     7, 3,     6, 4,     7, 4}
@@ -5207,7 +5229,7 @@ void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, doub
 
 
   xi=matrix(1,npar,1,npar);
-  for (i=1;i<=npar;i++)
+  for (i=1;i<=npar;i++)  /* Starting with canonical directions j=1,n xi[i=1,n][j] */
     for (j=1;j<=npar;j++)
       xi[i][j]=(i==j ? 1.0 : 0.0);
   printf("Powell\n");  fprintf(ficlog,"Powell\n");