[BACK]Return to imachprax.c CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/imachprax.c between versions 1.2 and 1.3

version 1.2, 2023/06/22 11:22:40 version 1.3, 2023/06/22 11:28:07
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.3  2023/06/22 11:28:07  brouard
     *** empty log message ***
   
   Revision 1.2  2023/06/22 11:22:40  brouard    Revision 1.2  2023/06/22 11:22:40  brouard
   Summary: with svd but not working yet    Summary: with svd but not working yet
   
Line 1281  Important routines Line 1284  Important routines
 /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */  /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
 /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */  /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
 /* #define FLATSUP  *//* Suppresses directions where likelihood is flat */  /* #define FLATSUP  *//* Suppresses directions where likelihood is flat */
   /* #define POWELLORIGINCONJUGATE  /\* Don't use conjugate but biggest decrease if valuable *\/ */
   
 #include <math.h>  #include <math.h>
 #include <stdio.h>  #include <stdio.h>
Line 2999  void powell(double p[], double **xi, int Line 3003  void powell(double p[], double **xi, int
       printf("%d",i);fflush(stdout); /* print direction (parameter) i */        printf("%d",i);fflush(stdout); /* print direction (parameter) i */
       fprintf(ficlog,"%d",i);fflush(ficlog);        fprintf(ficlog,"%d",i);fflush(ficlog);
 #ifdef LINMINORIGINAL  #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  #else
       linmin(p,xit,n,fret,func,&flat); /* Point p[n]. xit[n] has been loaded for direction i as input.*/        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 */                          flatdir[i]=flat; /* Function is vanishing in that direction i */
Line 3029  void powell(double p[], double **xi, int Line 3035  void powell(double p[], double **xi, int
       fprintf(ficlog,"\n");        fprintf(ficlog,"\n");
 #endif  #endif
     } /* end loop on each direction i */      } /* 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  */      /* 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++) {      for(j=1;j<=n;j++) {
       if(flatdir[j] >0){        if(flatdir[j] >0){
         printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);          printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
Line 3169  void powell(double p[], double **xi, int Line 3174  void powell(double p[], double **xi, int
           }            }
         }          }
 #endif  #endif
   #ifdef POWELLORIGINCONJUGATE
         for (j=1;j<=n;j++) {           for (j=1;j<=n;j++) { 
           xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */            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 */            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  #ifdef LINMINORIGINAL
 #else  #else
         for (j=1, flatd=0;j<=n;j++) {          for (j=1, flatd=0;j<=n;j++) {
Line 3201  void powell(double p[], double **xi, int Line 3213  void powell(double p[], double **xi, int
 #endif  #endif
         printf("Gaining to use new 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);          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  #ifdef DEBUG
         printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);          printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
Line 4811  double funcone( double *x) Line 4832  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 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 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          * 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[]          * 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 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]          * 3 we want to compute =cotvar[mw[mi][i]][TvarVVA[ncovva]][i] at position TvarVVAind[ncovva]
Line 4825  double funcone( double *x) Line 4846  double funcone( double *x)
         *                   2, 3, 4, 6, 7,          *                   2, 3, 4, 6, 7,
         *                     6, 8, 9, 10, 11}          *                     6, 8, 9, 10, 11}
         * TvarFind[itv]                        0      0       0          * 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[ncovf]]=[1]=2 [2]=3 [4]=4
         * Tvar[TvarFind[itv]]                [0]=?      ?ncovv 1 à ncovvt]          * 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}          *   Not a fixed cotvar[mw][itv][i]     6       7      6  2      7, 2,     6, 3,     7, 3,     6, 4,     7, 4}
Line 5210  void mlikeli(FILE *ficres,double p[], in Line 5232  void mlikeli(FILE *ficres,double p[], in
   
   
   xi=matrix(1,npar,1,npar);    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++)      for (j=1;j<=npar;j++)
       xi[i][j]=(i==j ? 1.0 : 0.0);        xi[i][j]=(i==j ? 1.0 : 0.0);
   printf("Powell\n");  fprintf(ficlog,"Powell\n");    printf("Powell\n");  fprintf(ficlog,"Powell\n");
Removed from v.1.2  
changed lines
  Added in v.1.3

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>