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

Diff for /imach/src/imach.c between versions 1.352 and 1.356

version 1.352, 2023/04/29 10:46:21 version 1.356, 2023/05/23 12:08:43
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.356  2023/05/23 12:08:43  brouard
     Summary: 0.99r46
   
     * imach.c (Module): Fixed PROB_r
   
     Revision 1.355  2023/05/22 17:03:18  brouard
     Summary: 0.99r46
   
     * imach.c (Module): In the ILK....txt file, the number of columns
     before the covariates values is dependent of the number of states (16+nlstate): 0.99r46
   
     Revision 1.354  2023/05/21 05:05:17  brouard
     Summary: Temporary change for imachprax
   
     Revision 1.353  2023/05/08 18:48:22  brouard
     *** empty log message ***
   
   Revision 1.352  2023/04/29 10:46:21  brouard    Revision 1.352  2023/04/29 10:46:21  brouard
   *** empty log message ***    *** empty log message ***
   
Line 1609  int **nbcode, *Tvar; /**< model=V2 => Tv Line 1626  int **nbcode, *Tvar; /**< model=V2 => Tv
 /* Tage[cptcovage]=k            5               8         10 */ /* Position in the model of ith cov*age */  /* Tage[cptcovage]=k            5               8         10 */ /* Position in the model of ith cov*age */
 /* model="V2+V3+V4+V6+V7+V6*V2+V7*V2+V6*V3+V7*V3+V6*V4+V7*V4+age*V2+age*V3+age*V4+age*V6+age*V7+age*V6*V2+age*V6*V3+age*V7*V3+age*V6*V4+age*V7*V4\r"*/  /* model="V2+V3+V4+V6+V7+V6*V2+V7*V2+V6*V3+V7*V3+V6*V4+V7*V4+age*V2+age*V3+age*V4+age*V6+age*V7+age*V6*V2+age*V6*V3+age*V7*V3+age*V6*V4+age*V7*V4\r"*/
 /*  p Tvard[1][1]@21 = {6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0}*/  /*  p Tvard[1][1]@21 = {6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0}*/
 /*  p Tvard[2][1]@21 = {7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0 <repeats 11 times>}  /*  p Tvard[2][1]@21 = {7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0 <repeats 11 times>} */
 /* p Tvardk[1][1]@24 = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0}*/  /* p Tvardk[1][1]@24 = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0}*/
 /* p Tvardk[1][1]@22 = {0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0} */  /* p Tvardk[1][1]@22 = {0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0} */
 /* Tvard[1][1]@4={4,3,1,2}    V4*V3 V1*V2               */ /* Position in model of the ith prod without age */  /* Tvard[1][1]@4={4,3,1,2}    V4*V3 V1*V2               */ /* Position in model of the ith prod without age */
Line 4550  double funcone( double *x) Line 4567  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=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 4565  double funcone( double *x) Line 4582  double funcone( double *x)
         *                     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       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 4576  double funcone( double *x) Line 4594  double funcone( double *x)
         ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/          ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
         /* if(TvarFind[itv]==0){ /\* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv *\/ */          /* if(TvarFind[itv]==0){ /\* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv *\/ */
         if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */          if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */
             /* printf("DEBUG ncovv=%d, Varying TvarVV[ncovv]=%d\n",ncovv, TvarVV[ncovv]); */
           cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i];  /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */             cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i];  /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */ 
             /* printf("DEBUG Varying cov[ioffset+ipos=%d]=%g \n",ioffset+ipos,cotvarv); */
         }else{ /* fixed covariate */          }else{ /* fixed covariate */
           /* cotvarv=covar[Tvar[TvarFind[itv]]][i];  /\* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model *\/ */            /* cotvarv=covar[Tvar[TvarFind[itv]]][i];  /\* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model *\/ */
             /* printf("DEBUG ncovv=%d, Fixed TvarVV[ncovv]=%d\n",ncovv, TvarVV[ncovv]); */
           cotvarv=covar[itv][i];  /* Good: In V6*V3, 3 is fixed at position of the data */            cotvarv=covar[itv][i];  /* Good: In V6*V3, 3 is fixed at position of the data */
             /* printf("DEBUG Fixed cov[ioffset+ipos=%d]=%g \n",ioffset+ipos,cotvarv); */
         }          }
         if(ipos!=iposold){ /* Not a product or first of a product */          if(ipos!=iposold){ /* Not a product or first of a product */
           cotvarvold=cotvarv;            cotvarvold=cotvarv;
Line 4588  double funcone( double *x) Line 4610  double funcone( double *x)
         }          }
         iposold=ipos;          iposold=ipos;
         cov[ioffset+ipos]=cotvarv;          cov[ioffset+ipos]=cotvarv;
           /* printf("DEBUG Product cov[ioffset+ipos=%d] \n",ioffset+ipos); */
         /* For products */          /* For products */
       }        }
       /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates single *\/ */        /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates single *\/ */
Line 4930  void mlikeli(FILE *ficres,double p[], in Line 4953  void mlikeli(FILE *ficres,double p[], in
   double fret;    double fret;
   double fretone; /* Only one call to likelihood */    double fretone; /* Only one call to likelihood */
   /*  char filerespow[FILENAMELENGTH];*/    /*  char filerespow[FILENAMELENGTH];*/
     
     double * p1; /* Shifted parameters from 0 instead of 1 */
 #ifdef NLOPT  #ifdef NLOPT
   int creturn;    int creturn;
   nlopt_opt opt;    nlopt_opt opt;
   /* double lb[9] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL }; /\* lower bounds *\/ */    /* double lb[9] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL }; /\* lower bounds *\/ */
   double *lb;    double *lb;
   double minf; /* the minimum objective value, upon return */    double minf; /* the minimum objective value, upon return */
   double * p1; /* Shifted parameters from 0 instead of 1 */  
   myfunc_data dinst, *d = &dinst;    myfunc_data dinst, *d = &dinst;
 #endif  #endif
   
Line 5272  double hessij( double x[], double **hess Line 5296  double hessij( double x[], double **hess
       kmax=kmax+10;        kmax=kmax+10;
     }      }
     if(kmax >=10 || firstime ==1){      if(kmax >=10 || firstime ==1){
         /* What are the thetai and thetaj? thetai/ncovmodel thetai=(thetai-thetai%ncovmodel)/ncovmodel +thetai%ncovmodel=(line,pos)  */
       printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);        printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
       fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);        fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
       printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);        printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
Line 7609  void varprob(char optionfilefiname[], do Line 7634  void varprob(char optionfilefiname[], do
    double ***varpij;     double ***varpij;
   
    strcpy(fileresprob,"PROB_");      strcpy(fileresprob,"PROB_"); 
    strcat(fileresprob,fileres);     strcat(fileresprob,fileresu);
    if((ficresprob=fopen(fileresprob,"w"))==NULL) {     if((ficresprob=fopen(fileresprob,"w"))==NULL) {
      printf("Problem with resultfile: %s\n", fileresprob);       printf("Problem with resultfile: %s\n", fileresprob);
      fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);       fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
Line 8357  true period expectancies (those weighted Line 8382  true period expectancies (those weighted
 /******************* Gnuplot file **************/  /******************* Gnuplot file **************/
 void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double bage, double fage , int prevfcast, int prevbcast, char pathc[], double p[], int offyear, int offbyear){  void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double bage, double fage , int prevfcast, int prevbcast, char pathc[], double p[], int offyear, int offbyear){
   
   char dirfileres[132],optfileres[132];    char dirfileres[256],optfileres[256];
   char gplotcondition[132], gplotlabel[132];    char gplotcondition[256], gplotlabel[256];
   int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,kf=0,kvar=0,kk=0,ipos=0,iposold=0,ij=0, ijp=0, l=0;    int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,kf=0,kvar=0,kk=0,ipos=0,iposold=0,ij=0, ijp=0, l=0;
   int lv=0, vlv=0, kl=0;    int lv=0, vlv=0, kl=0;
   int ng=0;    int ng=0;
Line 8429  void printinggnuplot(char fileresu[], ch Line 8454  void printinggnuplot(char fileresu[], ch
     kvar=Tvar[TvarFind[kf]]; /* variable name */      kvar=Tvar[TvarFind[kf]]; /* variable name */
     /* k=18+Tvar[TvarFind[kf]];/\*offset because there are 18 columns in the ILK_ file but could be placed else where *\/ */      /* k=18+Tvar[TvarFind[kf]];/\*offset because there are 18 columns in the ILK_ file but could be placed else where *\/ */
     /* k=18+kf;/\*offset because there are 18 columns in the ILK_ file *\/ */      /* k=18+kf;/\*offset because there are 18 columns in the ILK_ file *\/ */
     k=19+kf;/*offset because there are 19 columns in the ILK_ file */      /* k=19+kf;/\*offset because there are 19 columns in the ILK_ file *\/ */
       k=16+nlstate+kf;/*offset because there are 19 columns in the ILK_ file, first cov Vn on col 21 with 4 living states */
     for (i=1; i<= nlstate ; i ++) {      for (i=1; i<= nlstate ; i ++) {
       fprintf(ficgp,"\nset out \"%s-p%dj-%d.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i,kvar);        fprintf(ficgp,"\nset out \"%s-p%dj-%d.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i,kvar);
       fprintf(ficgp,"unset log;\n# For each simple dummy covariate of the model \n plot  \"%s\"",subdirf(fileresilk));        fprintf(ficgp,"unset log;\n# For each simple dummy covariate of the model \n plot  \"%s\"",subdirf(fileresilk));
Line 11122  int decoderesult( char resultline[], int Line 11148  int decoderesult( char resultline[], int
   if (strlen(resultsav) >1){    if (strlen(resultsav) >1){
     j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' in this resultline */      j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' in this resultline */
   }    }
   if(j == 0){ /* Resultline but no = */    if(j == 0 && cptcovs== 0){ /* Resultline but no =  and no covariate in the model */
     TKresult[nres]=0; /* Combination for the nresult and the model */      TKresult[nres]=0; /* Combination for the nresult and the model */
     return (0);      return (0);
   }    }
   if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */    if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */
     printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model);      fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, 1+age+%s.\n",j, cptcovs, model);fflush(ficlog);
     fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model);      printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, 1+age+%s.\n",j, cptcovs, model);fflush(stdout);
     /* return 1;*/      if(j==0)
         return 1;
   }    }
   for(k=1; k<=j;k++){ /* Loop on any covariate of the RESULT LINE */    for(k=1; k<=j;k++){ /* Loop on any covariate of the RESULT LINE */
     if(nbocc(resultsav,'=') >1){      if(nbocc(resultsav,'=') >1){
Removed from v.1.352  
changed lines
  Added in v.1.356

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