Diff for /imach/src/imach.c between versions 1.339 and 1.345

version 1.339, 2022/09/09 17:55:22 version 1.345, 2022/09/16 13:40:11
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.345  2022/09/16 13:40:11  brouard
     Summary: Version 0.99r41
   
     * imach.c (Module): 0.99r41 Was an error when product of timevarying and fixed. Using FixedV[of name] now. Thank you  Feinuo
   
     Revision 1.344  2022/09/14 19:33:30  brouard
     Summary: version 0.99r40
   
     * imach.c (Module): Fixing names of variables in T_ (thanks to Feinuo)
   
     Revision 1.343  2022/09/14 14:22:16  brouard
     Summary: version 0.99r39
   
     * imach.c (Module): Version 0.99r39 with colored dummy covariates
     (fixed or time varying), using new last columns of
     ILK_parameter.txt file.
   
     Revision 1.342  2022/09/11 19:54:09  brouard
     Summary: 0.99r38
   
     * imach.c (Module): Adding timevarying products of any kinds,
     should work before shifting cotvar from ncovcol+nqv columns in
     order to have a correspondance between the column of cotvar and
     the id of column.
     (Module): Some cleaning and adding covariates in ILK.txt
   
     Revision 1.341  2022/09/11 07:58:42  brouard
     Summary: Version 0.99r38
   
     After adding change in cotvar.
   
     Revision 1.340  2022/09/11 07:53:11  brouard
     Summary: Version imach 0.99r37
   
     * imach.c (Module): Adding timevarying products of any kinds,
     should work before shifting cotvar from ncovcol+nqv columns in
     order to have a correspondance between the column of cotvar and
     the id of column.
   
   Revision 1.339  2022/09/09 17:55:22  brouard    Revision 1.339  2022/09/09 17:55:22  brouard
   Summary: version 0.99r37    Summary: version 0.99r37
   
Line 1272  typedef struct { Line 1311  typedef struct {
 #define MAXLINE 2048 /* Was 256 and 1024. Overflow with 312 with 2 states and 4 covariates. Should be ok */  #define MAXLINE 2048 /* Was 256 and 1024. Overflow with 312 with 2 states and 4 covariates. Should be ok */
   
 #define GNUPLOTPROGRAM "gnuplot"  #define GNUPLOTPROGRAM "gnuplot"
   #define GNUPLOTVERSION 5.1
   double gnuplotversion=GNUPLOTVERSION;
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/  /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
 #define FILENAMELENGTH 256  #define FILENAMELENGTH 256
   
Line 1317  char fullversion[]="$Revision$ $Date$"; Line 1358  char fullversion[]="$Revision$ $Date$";
 char strstart[80];  char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];  char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */  int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
   int debugILK=0; /* debugILK is set by a #d in a comment line */
 int nagesqr=0, nforce=0; /* nagesqr=1 if model is including age*age, number of forces */  int nagesqr=0, nforce=0; /* nagesqr=1 if model is including age*age, number of forces */
 /* Number of covariates model (1)=V2+V1+ V3*age+V2*V4 */  /* Number of covariates model (1)=V2+V1+ V3*age+V2*V4 */
 /* Model(2)  V1 + V2 + V3 + V8 + V7*V8 + V5*V6 + V8*age + V3*age + age*age */  /* Model(2)  V1 + V2 + V3 + V8 + V7*V8 + V5*V6 + V8*age + V3*age + age*age */
Line 1495  double  **covar; /**< covar[j,i], value Line 1537  double  **covar; /**< covar[j,i], value
                   * covar=matrix(0,NCOVMAX,1,n);                     * covar=matrix(0,NCOVMAX,1,n); 
                   * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */                    * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */
 double **coqvar; /* Fixed quantitative covariate nqv */  double **coqvar; /* Fixed quantitative covariate nqv */
 double ***cotvar; /* Time varying covariate ntv */  double ***cotvar; /* Time varying covariate start at ncovcol + nqv + (1 to ntv) */
 double ***cotqvar; /* Time varying quantitative covariate itqv */  double ***cotqvar; /* Time varying quantitative covariate itqv */
 double  idx;   double  idx; 
 int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */  int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
Line 1507  int **nbcode, *Tvar; /**< model=V2 => Tv Line 1549  int **nbcode, *Tvar; /**< model=V2 => Tv
        *  cptcovn number of covariates (not including constant and age or age*age) = number of plus sign + 1 = 10+1=11         *  cptcovn number of covariates (not including constant and age or age*age) = number of plus sign + 1 = 10+1=11
        * For time varying covariate, quanti or dummies         * For time varying covariate, quanti or dummies
        *       cotqvar[wav][iv(1 to nqtv)][i]= [1][12][i]=(V12) quanti         *       cotqvar[wav][iv(1 to nqtv)][i]= [1][12][i]=(V12) quanti
        *       cotvar[wav][ntv+iv][i]= [3+(1 to nqtv)][i]=(V12) quanti         *       cotvar[wav][ncovcol+nqv+ iv(1 to nqtv)][i]= [(1 to nqtv)][i]=(V12) quanti
        *       cotvar[wav][iv(1 to ntv)][i]= [1][1][i]=(V9) dummies at wav 1         *       cotvar[wav][iv(1 to ntv)][i]= [1][1][i]=(V9) dummies at wav 1
        *       cotvar[wav][iv(1 to ntv)][i]= [1][2][i]=(V10) dummies at wav 1         *       cotvar[wav][iv(1 to ntv)][i]= [1][2][i]=(V10) dummies at wav 1
        *       covar[Vk,i], value of the Vkth fixed covariate dummy or quanti for individual i:         *       covar[Vk,i], value of the Vkth fixed covariate dummy or quanti for individual i:
Line 1520  int **nbcode, *Tvar; /**< model=V2 => Tv Line 1562  int **nbcode, *Tvar; /**< model=V2 => Tv
   # States 1=Coresidence, 2 Living alone, 3 Institution    # States 1=Coresidence, 2 Living alone, 3 Institution
   # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi    # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi
 */  */
 /*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */  /*           V5+V4+ V3+V4*V3 +V5*age+V2 +V1*V2+V1*age+V1 */
 /*    k        1  2   3   4     5    6    7     8    9 */  /*    kmodel  1  2   3   4     5    6    7     8    9 */
 /*Typevar[k]=  0  0   0   2     1    0    2     1    0 *//*0 for simple covariate (dummy, quantitative,*/  /*Typevar[k]=  0  0   0   2     1    0    2     1    0 *//*0 for simple covariate (dummy, quantitative,*/
                                                          /* fixed or varying), 1 for age product, 2 for*/                                                           /* fixed or varying), 1 for age product, 2 for*/
                                                          /* product */                                                           /* product */
Line 3880  double func( double *x) Line 3922  double func( double *x)
   int ioffset=0;    int ioffset=0;
   int ipos=0,iposold=0,ncovv=0;    int ipos=0,iposold=0,ncovv=0;
   
     double cotvarv, cotvarvold;
   double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];    double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
   double **out;    double **out;
   double lli; /* Individual log likelihood */    double lli; /* Individual log likelihood */
Line 3911  double func( double *x) Line 3954  double func( double *x)
       */        */
       ioffset=2+nagesqr ;        ioffset=2+nagesqr ;
    /* Fixed */     /* Fixed */
       for (kf=1; kf<=ncovf;kf++){ /* For each fixed covariate dummu or quant or prod */        for (kf=1; kf<=ncovf;kf++){ /* For each fixed covariate dummy or quant or prod */
         /* # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi */          /* # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi */
         /*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */          /*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
         /*  TvarF[1]=Tvar[6]=2,  TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1  ID of fixed covariates or product V2, V1*V2, V1 */          /*  TvarF[1]=Tvar[6]=2,  TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1  ID of fixed covariates or product V2, V1*V2, V1 */
Line 3927  double func( double *x) Line 3970  double func( double *x)
          mw[mi][i] is real wave of the mi th effectve wave */           mw[mi][i] is real wave of the mi th effectve wave */
       /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i];        /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i];
          s2=s[mw[mi+1][i]][i];           s2=s[mw[mi+1][i]][i];
          And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i]           And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i] because now is moved after nvocol+nqv 
          But if the variable is not in the model TTvar[iv] is the real variable effective in the model:           But if the variable is not in the model TTvar[iv] is the real variable effective in the model:
          meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]           meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]
       */        */
Line 3937  double func( double *x) Line 3980  double func( double *x)
         /*   /\* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; but where is the crossproduct? *\/ */          /*   /\* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; but where is the crossproduct? *\/ */
         /*   cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; */          /*   cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; */
         /* } */          /* } */
         for(ncovv=1, ipos=0; ncovv <= ncovvt ; ncovv++){ /* Varying  covariates (single and product but no age )*/          for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Varying  covariates (single and product but no age )*/
           itv=TvarVV[ncovv]; /*  TvarVV={3, 1, 3}  */            itv=TvarVV[ncovv]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate */
           ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] */            ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
             if(FixedV[itv]!=0){ /* Not a fixed covariate */
               cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i];  /* cotvar[wav][ncovcol+nqv+iv][i] */
             }else{ /* fixed covariate */
               cotvarv=covar[itv][i];  /* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model */
             }
           if(ipos!=iposold){ /* Not a product or first of a product */            if(ipos!=iposold){ /* Not a product or first of a product */
             /* TvarFind={1,0,0,0}  */              cotvarvold=cotvarv;
             if(TvarFind[itv]==0){            }else{ /* A second product */
             cov[ioffset+ipos]= cotvar[mw[mi][i]][ncovv][i];  /* Should be covar if fixed covar[Tvar[TvarFind[itv]]][i]*/              cotvarv=cotvarv*cotvarvold;
             }else{  
               cov[ioffset+ipos]=covar[Tvar[TvarFind[itv]]][i];  
             }  
           }else{  
             if(TvarFind[itv]==0){  
               cov[ioffset+ipos]*= cotvar[mw[mi][i]][ncovv][i];  /* Should be covar if fixed covar[Tvar[TvarFind[itv]]][i]*/  
             }else{  
               cov[ioffset+ipos]*=covar[Tvar[TvarFind[itv]]][i];  
             }  
           }            }
           iposold=ipos;            iposold=ipos;
           /* For products */            cov[ioffset+ipos]=cotvarv;
         }          }
         /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */          /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */
         /*   iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */          /*   iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */
Line 3988  double func( double *x) Line 4027  double func( double *x)
             if(!FixedV[Tvar[Tage[kk]]])              if(!FixedV[Tvar[Tage[kk]]])
               cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */                cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
             else              else
               cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact;                cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ 
           }            }
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,            out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                        1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));                         1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
Line 4084  double func( double *x) Line 4123  double func( double *x)
         }           } 
         /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/          /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
         /*if(lli ==000.0)*/          /*if(lli ==000.0)*/
         /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */          /* printf("num[i], i=%d, bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
         ipmx +=1;          ipmx +=1;
         sw += weight[i];          sw += weight[i];
         ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;          ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
Line 4101  double func( double *x) Line 4140  double func( double *x)
         cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];          cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];
       for(mi=1; mi<= wav[i]-1; mi++){        for(mi=1; mi<= wav[i]-1; mi++){
         for(k=1; k <= ncovv ; k++){          for(k=1; k <= ncovv ; k++){
           cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];            cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ 
         }          }
         for (ii=1;ii<=nlstate+ndeath;ii++)          for (ii=1;ii<=nlstate+ndeath;ii++)
           for (j=1;j<=nlstate+ndeath;j++){            for (j=1;j<=nlstate+ndeath;j++){
Line 4148  double func( double *x) Line 4187  double func( double *x)
           if(nagesqr==1)            if(nagesqr==1)
             cov[3]= agexact*agexact;              cov[3]= agexact*agexact;
           for (kk=1; kk<=cptcovage;kk++) {            for (kk=1; kk<=cptcovage;kk++) {
             cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;              if(!FixedV[Tvar[Tage[kk]]])
                 cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
               else
                 cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ 
           }            }
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,            out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                        1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));                         1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
Line 4204  double func( double *x) Line 4246  double func( double *x)
         ipmx +=1;          ipmx +=1;
         sw += weight[i];          sw += weight[i];
         ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;          ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
 /*      printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */          /* printf("num[i]=%09ld, i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
       } /* end of wave */        } /* end of wave */
     } /* end of individual */      } /* end of individual */
   }else{  /* ml=5 no inter-extrapolation no jackson =0.8a */    }else{  /* ml=5 no inter-extrapolation no jackson =0.8a */
Line 4223  double func( double *x) Line 4265  double func( double *x)
           if(nagesqr==1)            if(nagesqr==1)
             cov[3]= agexact*agexact;              cov[3]= agexact*agexact;
           for (kk=1; kk<=cptcovage;kk++) {            for (kk=1; kk<=cptcovage;kk++) {
             cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;              if(!FixedV[Tvar[Tage[kk]]])
                 cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
               else
                 cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ 
           }            }
                   
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,            out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
Line 4256  double funcone( double *x) Line 4301  double funcone( double *x)
   int ioffset=0;    int ioffset=0;
   int ipos=0,iposold=0,ncovv=0;    int ipos=0,iposold=0,ncovv=0;
   
     double cotvarv, cotvarvold;
   double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];    double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
   double **out;    double **out;
   double lli; /* Individual log likelihood */    double lli; /* Individual log likelihood */
Line 4310  double funcone( double *x) Line 4356  double funcone( double *x)
          mw[mi][i] is real wave of the mi th effectve wave */           mw[mi][i] is real wave of the mi th effectve wave */
       /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i];        /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i];
          s2=s[mw[mi+1][i]][i];           s2=s[mw[mi+1][i]][i];
          And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i]           And the iv th varying covariate in the DATA is the cotvar[mw[mi+1][i]][ncovcol+nqv+iv][i]
          But if the variable is not in the model TTvar[iv] is the real variable effective in the model:  
          meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]  
       */        */
     /* This part may be useless now because everythin should be in covar */      /* This part may be useless now because everythin should be in covar */
     /* for (k=1; k<=nqfveff;k++){ /\* Simple and product fixed Quantitative covariates without age* products *\/ */      /* for (k=1; k<=nqfveff;k++){ /\* Simple and product fixed Quantitative covariates without age* products *\/ */
Line 4337  double funcone( double *x) Line 4381  double funcone( double *x)
       /* We need the position of the time varying or product in the model */        /* We need the position of the time varying or product in the model */
       /* TvarVVind={2,5,5}, for V3 at position 2 and then the product V1*V3 is decomposed into V1 and V3 but at same position 5 */                     /* TvarVVind={2,5,5}, for V3 at position 2 and then the product V1*V3 is decomposed into V1 and V3 but at same position 5 */             
       /* TvarVV gives the variable name */        /* TvarVV gives the variable name */
       for(ncovv=1, ipos=0; ncovv <= ncovvt ; ncovv++){ /* Varying  covariates (single and product but no age )*/        /* Other example V1 + V3 + V5 + age*V1  + age*V3 + age*V5 + V1*V3  + V3*V5  + V1*V5 
         itv=TvarVV[ncovv]; /*  TvarVV={3, 1, 3}  */        *      k=         1   2     3     4         5        6        7       8        9
         ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] */        *  varying            1     2                                 3       4        5
         *  ncovv              1     2                                3 4     5 6      7 8
         * TvarVV[ncovv]      V3     5                                1 3     3 5      1 5
         * TvarVVind           2     3                                7 7     8 8      9 9
         * TvarFind[k]     1   0     0     0         0        0        0       0        0
         */
         /* Other model ncovcol=5 nqv=0 ntv=3 nqtv=0 nlstate=3
          /* V2 V3 V4 are fixed V6 V7 are timevarying so V8 and V5 are not in the model and product column will start at 9 Tvar[4]=6
           * FixedV[ncovcol+qv+ntv+nqtv]       V5
           *             V1  V2     V3    V4   V5 V6     V7  V8
           *             0   0      0      0    0  1      1   1 
           * model=          V2  +  V3  +  V4  +  V6  +  V7  +  V6*V2  +  V7*V2  +  V6*V3  +  V7*V3  +  V6*V4  +  V7*V4
           * kmodel           1     2      3      4      5        6         7         8         9        10        11
           * ncovf            1     2      3
           * ncovvt=14                            1      2       3 4       5 6       7 8       9 10     11 12     13 14
           * TvarVV[1]@14 = itv                   {6,     7,     6, 2,     7, 2,     6, 3,     7, 3,     6, 4,     7, 4}
           * TvarVVind[1]@14=                    {4,     5,      6, 6,     7, 7,     8, 8,      9, 9,   10, 10,   11, 11}
           * TvarFind[1]@14= {1,    2,     3,     0 <repeats 12 times>}
           * Tvar[1]@20=     {2,    3,     4,    6,      7,      9,      10,        11,       12,      13,       14}
           * TvarFind[itv]                        0      0       0
           * FixedV[itv]                          1      1       1  0      1 0       1 0       1 0       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}
           *   fixed covar[itv]                  [6]     [7]    [6][2]                            
           */
   
         for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Varying  covariates (single and product but no age) including individual from products */
           itv=TvarVV[ncovv]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate, exploding product  */
           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(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */
             cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i];  /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */ 
           }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[itv][i];  /* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model */
           }
         if(ipos!=iposold){ /* Not a product or first of a product */          if(ipos!=iposold){ /* Not a product or first of a product */
           /* TvarFind={1,0,0,0}  */            cotvarvold=cotvarv;
           if(TvarFind[itv]==0){          }else{ /* A second product */
             cov[ioffset+ipos]= cotvar[mw[mi][i]][ncovv][i];  /* Should be covar if fixed covar[Tvar[TvarFind[itv]]][i]*/            cotvarv=cotvarv*cotvarvold;
           }else{  
             cov[ioffset+ipos]=covar[Tvar[TvarFind[itv]]][i];  
           }  
         }else{  
           if(TvarFind[itv]==0){  
             cov[ioffset+ipos]*= cotvar[mw[mi][i]][ncovv][i];  /* Should be covar if fixed covar[Tvar[TvarFind[itv]]][i]*/  
           }else{  
             cov[ioffset+ipos]*=covar[Tvar[TvarFind[itv]]][i];  
           }  
         }          }
         iposold=ipos;          iposold=ipos;
           cov[ioffset+ipos]=cotvarv;
         /* 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 4397  double funcone( double *x) Line 4469  double funcone( double *x)
           if(!FixedV[Tvar[Tage[kk]]])            if(!FixedV[Tvar[Tage[kk]]])
             cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;              cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
           else            else
             cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact;              cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ 
         }          }
         /* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */          /* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */
         /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */          /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
Line 4452  double funcone( double *x) Line 4524  double funcone( double *x)
       ipmx +=1;        ipmx +=1;
       sw += weight[i];        sw += weight[i];
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;        ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
       /* printf("Funcone i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */        /* Printing covariates values for each contribution for checking */
         /* printf("num[i]=%09ld, i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
       if(globpr){        if(globpr){
         fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\          fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
  %11.6f %11.6f %11.6f ", \   %11.6f %11.6f %11.6f ", \
                 num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,                  num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
                 2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2]));                  2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2]));
  /*     printf("%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ */          /*      printf("%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ */
  /* %11.6f %11.6f %11.6f ", \ */          /* %11.6f %11.6f %11.6f ", \ */
  /*             num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, */          /*              num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, */
  /*             2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */          /*              2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
         for(k=1,llt=0.,l=0.; k<=nlstate; k++){          for(k=1,llt=0.,l=0.; k<=nlstate; k++){
           llt +=ll[k]*gipmx/gsw;            llt +=ll[k]*gipmx/gsw;
           fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);            fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
           /* printf(" %10.6f",-ll[k]*gipmx/gsw); */            /* printf(" %10.6f",-ll[k]*gipmx/gsw); */
         }          }
         fprintf(ficresilk," %10.6f\n", -llt);          fprintf(ficresilk," %10.6f ", -llt);
         /* printf(" %10.6f\n", -llt); */          /* printf(" %10.6f\n", -llt); */
       }          /* if(debugILK){ /\* debugILK is set by a #d in a comment line *\/ */
           /* fprintf(ficresilk,"%09ld ", num[i]); */ /* not necessary */
           for (kf=1; kf<=ncovf;kf++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */
             fprintf(ficresilk," %g",covar[Tvar[TvarFind[kf]]][i]);
           }
           for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Varying  covariates (single and product but no age) including individual from products */
             ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
             if(ipos!=iposold){ /* Not a product or first of a product */
               fprintf(ficresilk," %g",cov[ioffset+ipos]);
               /* printf(" %g",cov[ioffset+ipos]); */
             }else{
               fprintf(ficresilk,"*");
               /* printf("*"); */
             }
             iposold=ipos;
           }
           for (kk=1; kk<=cptcovage;kk++) {
             if(!FixedV[Tvar[Tage[kk]]]){
               fprintf(ficresilk," %g*age",covar[Tvar[Tage[kk]]][i]);
               /* printf(" %g*age",covar[Tvar[Tage[kk]]][i]); */
             }else{
               fprintf(ficresilk," %g*age",cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]);/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ 
               /* printf(" %g*age",cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]);/\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\/  */
             }
           }
           /* printf("\n"); */
           /* } /\*  End debugILK *\/ */
           fprintf(ficresilk,"\n");
         } /* End if globpr */
     } /* end of wave */      } /* end of wave */
   } /* end of individual */    } /* end of individual */
   for(k=1,l=0.; k<=nlstate; k++) l += ll[k];    for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
Line 4479  double funcone( double *x) Line 4580  double funcone( double *x)
     gipmx=ipmx;      gipmx=ipmx;
     gsw=sw;      gsw=sw;
   }    }
 return -l;    return -l;
 }  }
   
   
Line 4490  void likelione(FILE *ficres,double p[], Line 4591  void likelione(FILE *ficres,double p[],
      the selection of individuals/waves and       the selection of individuals/waves and
      to check the exact contribution to the likelihood.       to check the exact contribution to the likelihood.
      Plotting could be done.       Plotting could be done.
    */    */
   int k;    void pstamp(FILE *ficres);
     int k, kf, kk, kvar, ncovv, iposold, ipos;
   
   if(*globpri !=0){ /* Just counts and sums, no printings */    if(*globpri !=0){ /* Just counts and sums, no printings */
     strcpy(fileresilk,"ILK_");       strcpy(fileresilk,"ILK_"); 
Line 4500  void likelione(FILE *ficres,double p[], Line 4602  void likelione(FILE *ficres,double p[],
       printf("Problem with resultfile: %s\n", fileresilk);        printf("Problem with resultfile: %s\n", fileresilk);
       fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);        fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);
     }      }
       pstamp(ficresilk);fprintf(ficresilk,"# model=1+age+%s\n",model);
     fprintf(ficresilk, "#individual(line's_record) count ageb ageend s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");      fprintf(ficresilk, "#individual(line's_record) count ageb ageend s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
     fprintf(ficresilk, "#num_i ageb agend i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav ");      fprintf(ficresilk, "#num_i ageb agend i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav ");
     /*  i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */      /*  i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */
     for(k=1; k<=nlstate; k++)       for(k=1; k<=nlstate; k++) 
       fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);        fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);
     fprintf(ficresilk," -2*gipw/gsw*weight*ll(total)\n");      fprintf(ficresilk," -2*gipw/gsw*weight*ll(total) ");
   }  
       /* if(debugILK){ /\* debugILK is set by a #d in a comment line *\/ */
         for(kf=1;kf <= ncovf; kf++){
           fprintf(ficresilk,"V%d",Tvar[TvarFind[kf]]);
           /* printf("V%d",Tvar[TvarFind[kf]]); */
         }
         for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){
           ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate */
           if(ipos!=iposold){ /* Not a product or first of a product */
             /* printf(" %d",ipos); */
             fprintf(ficresilk," V%d",TvarVV[ncovv]);
           }else{
             /* printf("*"); */
             fprintf(ficresilk,"*");
           }
           iposold=ipos;
         }
         for (kk=1; kk<=cptcovage;kk++) {
           if(!FixedV[Tvar[Tage[kk]]]){
             /* printf(" %d*age(Fixed)",Tvar[Tage[kk]]); */
             fprintf(ficresilk," %d*age(Fixed)",Tvar[Tage[kk]]);
           }else{
             fprintf(ficresilk," %d*age(Varying)",Tvar[Tage[kk]]);/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ 
             /* printf(" %d*age(Varying)",Tvar[Tage[kk]]);/\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\/  */
           }
         }
       /* } /\* End if debugILK *\/ */
       /* printf("\n"); */
       fprintf(ficresilk,"\n");
     } /* End glogpri */
   
   *fretone=(*func)(p);    *fretone=(*func)(p);
   if(*globpri !=0){    if(*globpri !=0){
Line 4518  void likelione(FILE *ficres,double p[], Line 4650  void likelione(FILE *ficres,double p[],
     fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));      fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
     fprintf(fichtm,"\n<br>Equation of the model: <b>model=1+age+%s</b><br>\n",model);       fprintf(fichtm,"\n<br>Equation of the model: <b>model=1+age+%s</b><br>\n",model); 
               
     for (k=1; k<= nlstate ; k++) {  
       fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \  
 <img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k);  
     }  
     fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \      fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \
 <img src=\"%s-ori.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));  <img src=\"%s-ori.png\">\n",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
     fprintf(fichtm,"<br>- and by state of destination <a href=\"%s-dest.png\">%s-dest.png</a><br> \      fprintf(fichtm,"<br>- and by state of destination <a href=\"%s-dest.png\">%s-dest.png</a><br> \
 <img src=\"%s-dest.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));  <img src=\"%s-dest.png\">\n",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
       
       for (k=1; k<= nlstate ; k++) {
         fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br>\n \
   <img src=\"%s-p%dj.png\">\n",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k);
         for(kf=1; kf <= ncovf; kf++){ /* For each simple dummy covariate of the model */
           /* kvar=Tvar[TvarFind[kf]]; */ /* variable */
           fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j with colored covariate V%d. Same dot size of all points but with a different color for transitions with dummy variable V%d=1 at beginning of transition (keeping former color for V%d=0): <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
   <img src=\"%s-p%dj-%d.png\">",k,k,Tvar[TvarFind[kf]],Tvar[TvarFind[kf]],Tvar[TvarFind[kf]],subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,Tvar[TvarFind[kf]]);
         }
         for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Loop on the time varying extended covariates (with extension of Vn*Vm */
           ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate */
           kvar=TvarVV[ncovv]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate */
           /* printf("DebugILK fichtm ncovv=%d, kvar=TvarVV[ncovv]=V%d, ipos=TvarVVind[ncovv]=%d, Dummy[ipos]=%d, Typevar[ipos]=%d\n", ncovv,kvar,ipos,Dummy[ipos],Typevar[ipos]); */
           if(ipos!=iposold){ /* Not a product or first of a product */
             /* fprintf(ficresilk," V%d",TvarVV[ncovv]); */
             /* printf(" DebugILK fichtm ipos=%d != iposold=%d\n", ipos, iposold); */
             if(Dummy[ipos]==0 && Typevar[ipos]==0){ /* Only if dummy time varying: Dummy(0, 1=quant singor prod without age,2 dummy*age, 3quant*age) Typevar (0 single, 1=*age,2=Vn*vm)  */
               fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j with colored time varying dummy covariate V%d. Same dot size of all points but with a different color for transitions with dummy variable V%d=1 at beginning of transition (keeping former color for V%d=0): <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
   <img src=\"%s-p%dj-%d.png\">",k,k,kvar,kvar,kvar,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,kvar);
             } /* End only for dummies time varying (single?) */
           }else{ /* Useless product */
             /* printf("*"); */
             /* fprintf(ficresilk,"*"); */ 
           }
           iposold=ipos;
         } /* For each time varying covariate */
       } /* End loop on states */
   
   /*     if(debugILK){ */
   /*       for(kf=1; kf <= ncovf; kf++){ /\* For each simple dummy covariate of the model *\/ */
   /*      /\* kvar=Tvar[TvarFind[kf]]; *\/ /\* variable *\/ */
   /*      for (k=1; k<= nlstate ; k++) { */
   /*        fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j with colored covariate V%. Same dot size of all points but with a different color for transitions with dummy variable V%d=1 at beginning of transition (keeping former color for V%d=0): <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \ */
   /* <img src=\"%s-p%dj-%d.png\">",k,k,Tvar[TvarFind[kf]],Tvar[TvarFind[kf]],Tvar[TvarFind[kf]],subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,Tvar[TvarFind[kf]]); */
   /*      } */
   /*       } */
   /*       for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /\* Loop on the time varying extended covariates (with extension of Vn*Vm *\/ */
   /*      ipos=TvarVVind[ncovv]; /\* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate *\/ */
   /*      kvar=TvarVV[ncovv]; /\*  TvarVV={3, 1, 3} gives the name of each varying covariate *\/ */
   /*      /\* printf("DebugILK fichtm ncovv=%d, kvar=TvarVV[ncovv]=V%d, ipos=TvarVVind[ncovv]=%d, Dummy[ipos]=%d, Typevar[ipos]=%d\n", ncovv,kvar,ipos,Dummy[ipos],Typevar[ipos]); *\/ */
   /*      if(ipos!=iposold){ /\* Not a product or first of a product *\/ */
   /*        /\* fprintf(ficresilk," V%d",TvarVV[ncovv]); *\/ */
   /*        /\* printf(" DebugILK fichtm ipos=%d != iposold=%d\n", ipos, iposold); *\/ */
   /*        if(Dummy[ipos]==0 && Typevar[ipos]==0){ /\* Only if dummy time varying: Dummy(0, 1=quant singor prod without age,2 dummy*age, 3quant*age) Typevar (0 single, 1=*age,2=Vn*vm)  *\/ */
   /*          for (k=1; k<= nlstate ; k++) { */
   /*            fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Same dot size of all points but with a different color for transitions with dummy variable V%d=1 at beginning of transition (keeping former color for V%d=0): <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \ */
   /* <img src=\"%s-p%dj-%d.png\">",k,k,kvar,kvar,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,kvar); */
   /*          } /\* End state *\/ */
   /*        } /\* End only for dummies time varying (single?) *\/ */
   /*      }else{ /\* Useless product *\/ */
   /*        /\* printf("*"); *\/ */
   /*        /\* fprintf(ficresilk,"*"); *\/  */
   /*      } */
   /*      iposold=ipos; */
   /*       } /\* For each time varying covariate *\/ */
   /*     }/\* End debugILK *\/ */
     fflush(fichtm);      fflush(fichtm);
   }    }/* End globpri */
   return;    return;
 }  }
   
Line 5249  Title=%s <br>Datafile=%s Firstpass=%d La Line 5433  Title=%s <br>Datafile=%s Firstpass=%d La
               if(anyvaryingduminmodel==1){ /* Some are varying covariates */                if(anyvaryingduminmodel==1){ /* Some are varying covariates */
                 for (z1=1; z1<=cptcoveff; z1++) {                  for (z1=1; z1<=cptcoveff; z1++) {
                   if( Fixed[Tmodelind[z1]]==1){                    if( Fixed[Tmodelind[z1]]==1){
                     iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;                      /* iv= Tvar[Tmodelind[z1]]-ncovcol-nqv; /\* Good *\/ */
                       iv= Tvar[Tmodelind[z1]]; /* Good *//* because cotvar starts now at first at ncovcol+nqv+ntv */ 
                     if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality. If covariate's                       if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality. If covariate's 
                                                                                       value is -1, we don't select. It differs from the                                                                                         value is -1, we don't select. It differs from the 
                                                                                       constant and age model which counts them. */                                                                                        constant and age model which counts them. */
Line 5330  Title=%s <br>Datafile=%s Firstpass=%d La Line 5515  Title=%s <br>Datafile=%s Firstpass=%d La
         fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");           fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); 
         fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");           fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); 
         fprintf(ficlog, "\n#********** Variable ");           fprintf(ficlog, "\n#********** Variable "); 
         for (z1=1; z1<=cptcovs; z1++){          for (z1=1; z1<=cptcoveff; z1++){
           if(!FixedV[Tvaraff[z1]]){            if(!FixedV[Tvaraff[z1]]){
             printf( "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);              printf( "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
             fprintf(ficresp, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);              fprintf(ficresp, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
Line 5786  void prevalence(double ***probs, double Line 5971  void prevalence(double ***probs, double
         /* Tvar[Tmodelind[z1]] is the n of Vn; n-ncovcol-nqv is the first time varying covariate or iv */          /* Tvar[Tmodelind[z1]] is the n of Vn; n-ncovcol-nqv is the first time varying covariate or iv */
         for (z1=1; z1<=cptcoveff; z1++){          for (z1=1; z1<=cptcoveff; z1++){
           if( Fixed[Tmodelind[z1]]==1){            if( Fixed[Tmodelind[z1]]==1){
             iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;              iv= Tvar[Tmodelind[z1]];/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ 
             if (cotvar[m][iv][i]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality */              if (cotvar[m][iv][i]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality */
               bool=0;                bool=0;
           }else if( Fixed[Tmodelind[z1]]== 0)  /* fixed */            }else if( Fixed[Tmodelind[z1]]== 0)  /* fixed */
Line 6095  void  concatwav(int wav[], int **dh, int Line 6280  void  concatwav(int wav[], int **dh, int
    /* Loop on covariates without age and products and no quantitative variable */     /* Loop on covariates without age and products and no quantitative variable */
    for (k=1; k<=cptcovt; k++) { /* cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */     for (k=1; k<=cptcovt; k++) { /* cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */
      for (j=-1; (j < maxncov); j++) Ndum[j]=0;       for (j=-1; (j < maxncov); j++) Ndum[j]=0;
      printf("Testing k=%d, cptcovt=%d\n",k, cptcovt);       /* printf("Testing k=%d, cptcovt=%d\n",k, cptcovt); */
      if(Dummy[k]==0 && Typevar[k] !=1 && Typevar[k] != 2){ /* Dummy covariate and not age product nor fixed product */        if(Dummy[k]==0 && Typevar[k] !=1 && Typevar[k] != 2){ /* Dummy covariate and not age product nor fixed product */ 
        switch(Fixed[k]) {         switch(Fixed[k]) {
        case 0: /* Testing on fixed dummy covariate, simple or product of fixed */         case 0: /* Testing on fixed dummy covariate, simple or product of fixed */
Line 6193  void  concatwav(int wav[], int **dh, int Line 6378  void  concatwav(int wav[], int **dh, int
          break;           break;
        } /* end switch */         } /* end switch */
      } /* end dummy test */       } /* end dummy test */
      if(Dummy[k]==1 && Typevar[k] !=1){ /* Quantitative covariate and not age product */        if(Dummy[k]==1 && Typevar[k] !=1 && Fixed ==0){ /* Fixed Quantitative covariate and not age product */ 
        for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the  modality of this covariate Vj*/         for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the  modality of this covariate Vj*/
          if(Tvar[k]<=0 || Tvar[k]>=NCOVMAX){           if(Tvar[k]<=0 || Tvar[k]>=NCOVMAX){
            printf("Error k=%d \n",k);             printf("Error k=%d \n",k);
Line 6659  void  concatwav(int wav[], int **dh, int Line 6844  void  concatwav(int wav[], int **dh, int
    /*   fprintf(ficresprobmorprev," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */     /*   fprintf(ficresprobmorprev," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */
    /* } */     /* } */
    for (j=1; j<= cptcovs; j++){ /* For each selected (single) quantitative value */ /* To be done*/     for (j=1; j<= cptcovs; j++){ /* For each selected (single) quantitative value */ /* To be done*/
        /* fprintf(ficresprobmorprev," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); */
      fprintf(ficresprobmorprev," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);       fprintf(ficresprobmorprev," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
    }     }
    /* for(j=1;j<=cptcoveff;j++)  */     /* for(j=1;j<=cptcoveff;j++)  */
Line 7293  To be simple, these graphs help to under Line 7479  To be simple, these graphs help to under
   
    for(nres=1;nres <=nresult; nres++){ /* For each resultline */     for(nres=1;nres <=nresult; nres++){ /* For each resultline */
    for(j1=1; j1<=tj;j1++){ /* For any combination of dummy covariates, fixed and varying */     for(j1=1; j1<=tj;j1++){ /* For any combination of dummy covariates, fixed and varying */
      printf("Varprob  TKresult[nres]=%d j1=%d, nres=%d, cptcovn=%d, cptcoveff=%d tj=%d cptcovs=%d\n",  TKresult[nres], j1, nres, cptcovn, cptcoveff, tj, cptcovs);       /* printf("Varprob  TKresult[nres]=%d j1=%d, nres=%d, cptcovn=%d, cptcoveff=%d tj=%d cptcovs=%d\n",  TKresult[nres], j1, nres, cptcovn, cptcoveff, tj, cptcovs); */
      if(tj != 1 && TKresult[nres]!= j1)       if(tj != 1 && TKresult[nres]!= j1)
        continue;         continue;
   
Line 7309  To be simple, these graphs help to under Line 7495  To be simple, these graphs help to under
   
        /* Including quantitative variables of the resultline to be done */         /* Including quantitative variables of the resultline to be done */
        for (z1=1; z1<=cptcovs; z1++){ /* Loop on each variable of this resultline  */         for (z1=1; z1<=cptcovs; z1++){ /* Loop on each variable of this resultline  */
          printf("Varprob modelresult[%d][%d]=%d model=1+age+%s \n",nres, z1, modelresult[nres][z1], model);           /* printf("Varprob modelresult[%d][%d]=%d model=1+age+%s \n",nres, z1, modelresult[nres][z1], model); */
          fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=1+age+%s \n",nres, z1, modelresult[nres][z1], model);           fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=1+age+%s \n",nres, z1, modelresult[nres][z1], model);
          /* fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=1+age+%s resultline[%d]=%s \n",nres, z1, modelresult[nres][z1], model, nres, resultline[nres]); */           /* fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=1+age+%s resultline[%d]=%s \n",nres, z1, modelresult[nres][z1], model, nres, resultline[nres]); */
          if(Dummy[modelresult[nres][z1]]==0){/* Dummy variable of the variable in position modelresult in the model corresponding to z1 in resultline  */           if(Dummy[modelresult[nres][z1]]==0){/* Dummy variable of the variable in position modelresult in the model corresponding to z1 in resultline  */
Line 7968  void printinggnuplot(char fileresu[], ch Line 8154  void printinggnuplot(char fileresu[], ch
   
   char dirfileres[132],optfileres[132];    char dirfileres[132],optfileres[132];
   char gplotcondition[132], gplotlabel[132];    char gplotcondition[132], gplotlabel[132];
   int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=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;
   int vpopbased;    int vpopbased;
Line 7994  void printinggnuplot(char fileresu[], ch Line 8180  void printinggnuplot(char fileresu[], ch
   fprintf(ficgp,"yoff=(%d > 2? 0:1);\n",nlstate);    fprintf(ficgp,"yoff=(%d > 2? 0:1);\n",nlstate);
   fprintf(ficgp,"\n#Peripheral arrows\nset for [i=1:%d] for [j=1:%d] arrow i*10+j from cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.95*(cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0) - cos(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta2:0)), -0.95*(sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) - sin(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d))+( i!=j?(i-j)/abs(i-j)*delta2:0)) ls (i < j? 1:2)\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);    fprintf(ficgp,"\n#Peripheral arrows\nset for [i=1:%d] for [j=1:%d] arrow i*10+j from cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.95*(cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0) - cos(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta2:0)), -0.95*(sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) - sin(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d))+( i!=j?(i-j)/abs(i-j)*delta2:0)) ls (i < j? 1:2)\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
   
   fprintf(ficgp,"\n#Centripete arrows (turning in other direction (1-i) instead of (i-1)) \nset for [i=1:%d] arrow (%d+1)*10+i from cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.80*(cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0)  ), -0.80*(sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) + yoff ) ls 4\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);    fprintf(ficgp,"\n#Centripete arrows (turning in other direction (1-i) instead of (i-1)) \nset for [i=1:%d] for [j=1:%d] arrow (%d+1)*10+i from cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.80*(cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0)  ), -0.80*(sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) + yoff ) ls 4\n",nlstate, nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
   fprintf(ficgp,"\n#show arrow\nunset label\n");    fprintf(ficgp,"\n#show arrow\nunset label\n");
   fprintf(ficgp,"\n#States labels, starting from 2 (2-i) instead of (1-i), was (i-1)\nset for [i=1:%d] label i sprintf(\"State %%d\",i) center at cos(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)), yoff+sin(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)) font \"helvetica, 16\" tc rgbcolor \"blue\"\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);    fprintf(ficgp,"\n#States labels, starting from 2 (2-i) instead of (1-i), was (i-1)\nset for [i=1:%d] label i sprintf(\"State %%d\",i) center at cos(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)), yoff+sin(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)) font \"helvetica, 16\" tc rgbcolor \"blue\"\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
   fprintf(ficgp,"\nset label %d+1 sprintf(\"State %%d\",%d+1) center at 0.,0.  font \"helvetica, 16\" tc rgbcolor \"red\"\n",nlstate,nlstate);    fprintf(ficgp,"\nset label %d+1 sprintf(\"State %%d\",%d+1) center at 0.,0.  font \"helvetica, 16\" tc rgbcolor \"red\"\n",nlstate,nlstate);
Line 8031  void printinggnuplot(char fileresu[], ch Line 8217  void printinggnuplot(char fileresu[], ch
   fprintf(ficgp,"\nset out;unset log\n");    fprintf(ficgp,"\nset out;unset log\n");
   /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */    /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
   
     /* Plot the probability implied in the likelihood by covariate value */
     fprintf(ficgp,"\nset ter pngcairo size 640, 480");
     /* if(debugILK==1){ */
     for(kf=1; kf <= ncovf; kf++){ /* For each simple dummy covariate of the model */
       kvar=Tvar[TvarFind[kf]]; /* variable */
       k=18+Tvar[TvarFind[kf]];/*offset because there are 18 columns in the ILK_ file */
       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,"unset log;\n# For each simple dummy covariate of the model \n plot  \"%s\"",subdirf(fileresilk));
         fprintf(ficgp,"  u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? 7 : 9):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt variable ps 0.4 lc variable \\\n",i,1,k,k,i,1,kvar);
         for (j=2; j<= nlstate+ndeath ; j ++) {
           fprintf(ficgp,",\\\n \"\" u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? 7 : 9):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt variable ps 0.4 lc variable ",i,j,k,k,i,j,kvar);
         }
         fprintf(ficgp,";\nset out; unset ylabel;\n"); 
       }
     } /* End of each covariate dummy */
     for(ncovv=1, iposold=0, kk=0; ncovv <= ncovvt ; ncovv++){
       /* Other example        V1 + V3 + V5 + age*V1  + age*V3 + age*V5 + V1*V3  + V3*V5  + V1*V5 
        *     kmodel       =     1   2     3     4         5        6        7       8        9
        *  varying                   1     2                                 3       4        5
        *  ncovv                     1     2                                3 4     5 6      7 8
        * TvarVV[ncovv]             V3     5                                1 3     3 5      1 5
        * TvarVVind[ncovv]=kmodel    2     3                                7 7     8 8      9 9
        * TvarFind[kmodel]       1   0     0     0         0        0        0       0        0
        * kdata     ncovcol=[V1 V2] nqv=0 ntv=[V3 V4] nqtv=V5
        * Dummy[kmodel]          0   0     1     2         2        3        1       1        1
        */
       ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate */
       kvar=TvarVV[ncovv]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate */
       /* printf("DebugILK ficgp ncovv=%d, kvar=TvarVV[ncovv]=%d, ipos=TvarVVind[ncovv]=%d, Dummy[ipos]=%d, Typevar[ipos]=%d\n", ncovv,kvar,ipos,Dummy[ipos],Typevar[ipos]); */
       if(ipos!=iposold){ /* Not a product or first of a product */
         /* printf(" %d",ipos); */
         /* fprintf(ficresilk," V%d",TvarVV[ncovv]); */
         /* printf(" DebugILK ficgp suite ipos=%d != iposold=%d\n", ipos, iposold); */
         kk++; /* Position of the ncovv column in ILK_ */
         k=18+ncovf+kk; /*offset because there are 18 columns in the ILK_ file plus ncovf fixed covariate */
         if(Dummy[ipos]==0 && Typevar[ipos]==0){ /* Only if dummy time varying: Dummy(0, 1=quant singor prod without age,2 dummy*age, 3quant*age) Typevar (0 single, 1=*age,2=Vn*vm)  */
           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,"unset log;\n# For each simple dummy covariate of the model \n plot  \"%s\"",subdirf(fileresilk));
   
             if(gnuplotversion >=5.2){ /* Former gnuplot versions do not have variable pointsize!! */
               /* printf("DebugILK gnuplotversion=%g >=5.2\n",gnuplotversion); */
               fprintf(ficgp,"  u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? 7 : 9):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt variable ps 0.4 lc variable \\\n",i,1,k,k,i,1,kvar);
               for (j=2; j<= nlstate+ndeath ; j ++) {
                 fprintf(ficgp,",\\\n \"\" u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? 7 : 9):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt variable ps 0.4 lc variable ",i,j,k,k,i,j,kvar);
               }
             }else{
               /* printf("DebugILK gnuplotversion=%g <5.2\n",gnuplotversion); */
               fprintf(ficgp,"  u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt 7 ps 0.4 lc variable \\\n",i,1,k,i,1,kvar);
               for (j=2; j<= nlstate+ndeath ; j ++) {
                 fprintf(ficgp,",\\\n \"\" u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt 7 ps 0.4 lc variable ",i,j,k,i,j,kvar);
               }
             }
             fprintf(ficgp,";\nset out; unset ylabel;\n"); 
           }
         }/* End if dummy varying */
       }else{ /*Product */
         /* printf("*"); */
         /* fprintf(ficresilk,"*"); */
       }
       iposold=ipos;
     } /* For each time varying covariate */
     /* } /\* debugILK==1 *\/ */
     /* unset log; plot  "rrtest1_sorted_4/ILK_rrtest1_sorted_4.txt" u  2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with points lc variable */                
     /* fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); */
     /* fprintf(ficgp,"\nreplot  \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); */
     fprintf(ficgp,"\nset out;unset log\n");
     /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
   
   
     
   strcpy(dirfileres,optionfilefiname);    strcpy(dirfileres,optionfilefiname);
   strcpy(optfileres,"vpl");    strcpy(optfileres,"vpl");
   /* 1eme*/    /* 1eme*/
Line 10317  int readdata(char datafile[], int firsto Line 10575  int readdata(char datafile[], int firsto
         if(strb[0]=='.') { /* Missing value */          if(strb[0]=='.') { /* Missing value */
           lval=-1;            lval=-1;
           cotqvar[j][iv][i]=-1; /* 0.0/0.0 */            cotqvar[j][iv][i]=-1; /* 0.0/0.0 */
           cotvar[j][ntv+iv][i]=-1; /* For performance reasons */            cotvar[j][ncovcol+nqv+ntv+iv][i]=-1; /* For performance reasons */
           if(isalpha(strb[1])) { /* .m or .d Really Missing value */            if(isalpha(strb[1])) { /* .m or .d Really Missing value */
             printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);              printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);
             fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);fflush(ficlog);              fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);fflush(ficlog);
Line 10337  int readdata(char datafile[], int firsto Line 10595  int readdata(char datafile[], int firsto
             return 1;              return 1;
           }            }
           cotqvar[j][iv][i]=dval;             cotqvar[j][iv][i]=dval; 
           cotvar[j][ntv+iv][i]=dval;             cotvar[j][ncovcol+nqv+ntv+iv][i]=dval; /* because cotvar starts now at first ntv */ 
         }          }
         strcpy(line,stra);          strcpy(line,stra);
       }/* end loop ntqv */        }/* end loop ntqv */
Line 10377  int readdata(char datafile[], int firsto Line 10635  int readdata(char datafile[], int firsto
  Exiting.\n",lval,linei, i,line,iv,j);fflush(ficlog);   Exiting.\n",lval,linei, i,line,iv,j);fflush(ficlog);
           return 1;            return 1;
         }          }
         cotvar[j][iv][i]=(double)(lval);          cotvar[j][ncovcol+nqv+iv][i]=(double)(lval);
         strcpy(line,stra);          strcpy(line,stra);
       }/* end loop ntv */        }/* end loop ntv */
               
Line 10590  int decoderesult( char resultline[], int Line 10848  int decoderesult( char resultline[], int
   printf("decoderesult:%s\n",resultline);    printf("decoderesult:%s\n",resultline);
   
   strcpy(resultsav,resultline);    strcpy(resultsav,resultline);
   printf("Decoderesult resultsav=\"%s\" resultline=\"%s\"\n", resultsav, resultline);    /* printf("Decoderesult resultsav=\"%s\" resultline=\"%s\"\n", resultsav, resultline); */
   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 */
   }    }
Line 10650  int decoderesult( char resultline[], int Line 10908  int decoderesult( char resultline[], int
         if(Tvar[k1]==Tvarsel[k2]) {/* Tvar is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5   */          if(Tvar[k1]==Tvarsel[k2]) {/* Tvar is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5   */
           modelresult[nres][k2]=k1;/* we found a Vn=1 corrresponding to Vn*age in the model modelresult[2]=1 modelresult[1]=2  modelresult[3]=3  modelresult[6]=4 modelresult[9]=5 */            modelresult[nres][k2]=k1;/* we found a Vn=1 corrresponding to Vn*age in the model modelresult[2]=1 modelresult[1]=2  modelresult[3]=3  modelresult[6]=4 modelresult[9]=5 */
           resultmodel[nres][k1]=k2; /* Added here */            resultmodel[nres][k1]=k2; /* Added here */
           printf("Decoderesult first modelresult[k2=%d]=%d (k1) V%d*AGE\n",k2,k1,Tvar[k1]);            /* printf("Decoderesult first modelresult[k2=%d]=%d (k1) V%d*AGE\n",k2,k1,Tvar[k1]); */
           match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */            match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */
           break;            break;
         }          }
Line 10663  int decoderesult( char resultline[], int Line 10921  int decoderesult( char resultline[], int
     }else if(Typevar[k1]==2){ /* Product No age We want to get the position in the resultline of the product in the model line*/      }else if(Typevar[k1]==2){ /* Product No age We want to get the position in the resultline of the product in the model line*/
       /* resultmodel[nres][of such a Vn * Vm product k1] is not unique, so can't exist, we feed Tvard[k1][1] and [2] */         /* resultmodel[nres][of such a Vn * Vm product k1] is not unique, so can't exist, we feed Tvard[k1][1] and [2] */ 
       match=0;        match=0;
       printf("Decoderesult very first Product Tvardk[k1=%d][1]=%d Tvardk[k1=%d][2]=%d V%d * V%d\n",k1,Tvardk[k1][1],k1,Tvardk[k1][2],Tvardk[k1][1],Tvardk[k1][2]);        /* printf("Decoderesult very first Product Tvardk[k1=%d][1]=%d Tvardk[k1=%d][2]=%d V%d * V%d\n",k1,Tvardk[k1][1],k1,Tvardk[k1][2],Tvardk[k1][1],Tvardk[k1][2]); */
       for(k2=1; k2 <=j;k2++){/* Loop on resultline. In result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */        for(k2=1; k2 <=j;k2++){/* Loop on resultline. In result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
         if(Tvardk[k1][1]==Tvarsel[k2]) {/* Tvardk is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5   */          if(Tvardk[k1][1]==Tvarsel[k2]) {/* Tvardk is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5   */
           /* modelresult[k2]=k1; */            /* modelresult[k2]=k1; */
           printf("Decoderesult first Product modelresult[k2=%d]=%d (k1) V%d * \n",k2,k1,Tvarsel[k2]);            /* printf("Decoderesult first Product modelresult[k2=%d]=%d (k1) V%d * \n",k2,k1,Tvarsel[k2]); */
           match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */            match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */
         }          }
       }        }
Line 10680  int decoderesult( char resultline[], int Line 10938  int decoderesult( char resultline[], int
       for(k2=1; k2 <=j;k2++){/* Loop on resultline. In result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */        for(k2=1; k2 <=j;k2++){/* Loop on resultline. In result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
         if(Tvardk[k1][2]==Tvarsel[k2]) {/* Tvardk is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5   */          if(Tvardk[k1][2]==Tvarsel[k2]) {/* Tvardk is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5   */
           /* modelresult[k2]=k1;*/            /* modelresult[k2]=k1;*/
           printf("Decoderesult second Product modelresult[k2=%d]=%d (k1) * V%d \n ",k2,k1,Tvarsel[k2]);            /* printf("Decoderesult second Product modelresult[k2=%d]=%d (k1) * V%d \n ",k2,k1,Tvarsel[k2]); */
           match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */            match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */
           break;            break;
         }          }
Line 10740  int decoderesult( char resultline[], int Line 10998  int decoderesult( char resultline[], int
     /* k counting number of combination of single dummies in the equation model */      /* k counting number of combination of single dummies in the equation model */
     /* k4 counting single dummies in the equation model */      /* k4 counting single dummies in the equation model */
     /* k4q counting single quantitatives in the equation model */      /* k4q counting single quantitatives in the equation model */
     if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Dummy and Single, k1 is sorting according to MODEL, but k3 to resultline */      if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Dummy and Single, fixed or timevarying, k1 is sorting according to MODEL, but k3 to resultline */
        /* k4+1= (not always if quant in model) position in the resultline V(Tvarsel)=Tvalsel=Tresult[nres][pos](value); V(Tvresult[nres][pos] (variable): V(variable)=value) */         /* k4+1= (not always if quant in model) position in the resultline V(Tvarsel)=Tvalsel=Tresult[nres][pos](value); V(Tvresult[nres][pos] (variable): V(variable)=value) */
       /* modelresult[k3]=k1: k3th position in the result line corresponds to the k1 position in the model line (doesn't work with products)*/        /* modelresult[k3]=k1: k3th position in the result line corresponds to the k1 position in the model line (doesn't work with products)*/
       /* Value in the (current nres) resultline of the variable at the k1th position in the model equation resultmodel[nres][k1]= k3 */        /* Value in the (current nres) resultline of the variable at the k1th position in the model equation resultmodel[nres][k1]= k3 */
Line 10764  int decoderesult( char resultline[], int Line 11022  int decoderesult( char resultline[], int
       Tvresult[nres][k3]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */        Tvresult[nres][k3]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */
       Tinvresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */        Tinvresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */
       precov[nres][k1]=Tvalsel[k3]; /* Value from resultline of the variable at the k1 position in the model */        precov[nres][k1]=Tvalsel[k3]; /* Value from resultline of the variable at the k1 position in the model */
       printf("Decoderesult Dummy k=%d, k1=%d precov[nres=%d][k1=%d]=%.f V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k1, nres, k1,precov[nres][k1], k2, k3, (int)Tvalsel[k3], k4);        /* printf("Decoderesult Dummy k=%d, k1=%d precov[nres=%d][k1=%d]=%.f V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k1, nres, k1,precov[nres][k1], k2, k3, (int)Tvalsel[k3], k4); */
       k4++;;        k4++;;
     }else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Quantitative and single */      }else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Quantitative and single */
       /* Tqresult[nres][result_position]= value of the variable at the result_position in the nres resultline                                 */        /* Tqresult[nres][result_position]= value of the variable at the result_position in the nres resultline                                 */
Line 10782  int decoderesult( char resultline[], int Line 11040  int decoderesult( char resultline[], int
       Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */        Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */
       TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */        TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */
       precov[nres][k1]=Tvalsel[k3q];        precov[nres][k1]=Tvalsel[k3q];
       printf("Decoderesult Quantitative nres=%d,precov[nres=%d][k1=%d]=%.f V(k2q=V%d)= Tvalsel[%d]=%d, Tvarsel[%d]=%f\n",nres, nres, k1,precov[nres][k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]);        /* printf("Decoderesult Quantitative nres=%d,precov[nres=%d][k1=%d]=%.f V(k2q=V%d)= Tvalsel[%d]=%d, Tvarsel[%d]=%f\n",nres, nres, k1,precov[nres][k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]); */
       k4q++;;        k4q++;;
     }else if( Dummy[k1]==2 ){ /* For dummy with age product */      }else if( Dummy[k1]==2 ){ /* For dummy with age product */
       /* Tvar[k1]; */ /* Age variable */        /* Tvar[k1]; */ /* Age variable */
Line 10792  int decoderesult( char resultline[], int Line 11050  int decoderesult( char resultline[], int
       k2=(int)Tvarsel[k3]; /* nres=1 k1=2=>k3=1 Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 (V4); k1=3=>k3=2 Tvarsel[2]=3 (V3)*/        k2=(int)Tvarsel[k3]; /* nres=1 k1=2=>k3=1 Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 (V4); k1=3=>k3=2 Tvarsel[2]=3 (V3)*/
       TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* TinvDoQresult[nres][4]=1 */        TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* TinvDoQresult[nres][4]=1 */
       precov[nres][k1]=Tvalsel[k3];        precov[nres][k1]=Tvalsel[k3];
       printf("Decoderesult Dummy with age k=%d, k1=%d precov[nres=%d][k1=%d]=%.f Tvar[%d]=V%d k2=Tvarsel[%d]=%d Tvalsel[%d]=%d\n",k, k1, nres, k1,precov[nres][k1], k1, Tvar[k1], k3,(int)Tvarsel[k3], k3, (int)Tvalsel[k3]);        /* printf("Decoderesult Dummy with age k=%d, k1=%d precov[nres=%d][k1=%d]=%.f Tvar[%d]=V%d k2=Tvarsel[%d]=%d Tvalsel[%d]=%d\n",k, k1, nres, k1,precov[nres][k1], k1, Tvar[k1], k3,(int)Tvarsel[k3], k3, (int)Tvalsel[k3]); */
     }else if( Dummy[k1]==3 ){ /* For quant with age product */      }else if( Dummy[k1]==3 ){ /* For quant with age product */
       k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 25.1=k3q */        k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 25.1=k3q */
       k2q=(int)Tvarsel[k3q]; /*  Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */        k2q=(int)Tvarsel[k3q]; /*  Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */
       TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* TinvDoQresult[nres][5]=25.1 */        TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* TinvDoQresult[nres][5]=25.1 */
       precov[nres][k1]=Tvalsel[k3q];        precov[nres][k1]=Tvalsel[k3q];
       printf("Decoderesult Quantitative with age nres=%d, k1=%d, precov[nres=%d][k1=%d]=%f Tvar[%d]=V%d V(k2q=%d)= Tvarsel[%d]=%d, Tvalsel[%d]=%f\n",nres, k1, nres, k1,precov[nres][k1], k1,  Tvar[k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]);        /* printf("Decoderesult Quantitative with age nres=%d, k1=%d, precov[nres=%d][k1=%d]=%f Tvar[%d]=V%d V(k2q=%d)= Tvarsel[%d]=%d, Tvalsel[%d]=%f\n",nres, k1, nres, k1,precov[nres][k1], k1,  Tvar[k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]); */
     }else if(Typevar[k1]==2 ){ /* For product quant or dummy (not with age) */      }else if(Typevar[k1]==2 ){ /* For product quant or dummy (not with age) */
       precov[nres][k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]];              precov[nres][k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]];      
       printf("Decoderesult Quantitative or Dummy (not with age) nres=%d k1=%d precov[nres=%d][k1=%d]=%.f V%d(=%.f) * V%d(=%.f) \n",nres, k1, nres, k1,precov[nres][k1], Tvardk[k1][1], TinvDoQresult[nres][Tvardk[k1][1]], Tvardk[k1][2], TinvDoQresult[nres][Tvardk[k1][2]]);        /* printf("Decoderesult Quantitative or Dummy (not with age) nres=%d k1=%d precov[nres=%d][k1=%d]=%.f V%d(=%.f) * V%d(=%.f) \n",nres, k1, nres, k1,precov[nres][k1], Tvardk[k1][1], TinvDoQresult[nres][Tvardk[k1][1]], Tvardk[k1][2], TinvDoQresult[nres][Tvardk[k1][2]]); */
     }else{      }else{
       printf("Error Decoderesult probably a product  Dummy[%d]==%d && Typevar[%d]==%d\n", k1, Dummy[k1], k1, Typevar[k1]);        printf("Error Decoderesult probably a product  Dummy[%d]==%d && Typevar[%d]==%d\n", k1, Dummy[k1], k1, Typevar[k1]);
       fprintf(ficlog,"Error Decoderesult probably a product  Dummy[%d]==%d && Typevar[%d]==%d\n", k1, Dummy[k1], k1, Typevar[k1]);        fprintf(ficlog,"Error Decoderesult probably a product  Dummy[%d]==%d && Typevar[%d]==%d\n", k1, Dummy[k1], k1, Typevar[k1]);
Line 10890  int decodemodel( char model[], int lasto Line 11148  int decodemodel( char model[], int lasto
        * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8         * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8
        *   k=  1    2      3       4     5       6      7        8         *   k=  1    2      3       4     5       6      7        8
        *  cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8         *  cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8
        *  covar[k,i], value of kth covariate if not including age for individual i:         *  covar[k,i], are for fixed covariates, value of kth covariate if not including age for individual i:
        *       covar[1][i]= (V1), covar[4][i]=(V4), covar[8][i]=(V8)         *       covar[1][i]= (V1), covar[4][i]=(V4), covar[8][i]=(V8)
        *  Tvar[k] # of the kth covariate:  Tvar[1]=2  Tvar[2]=1 Tvar[4]=3 Tvar[8]=8         *  Tvar[k] # of the kth covariate:  Tvar[1]=2  Tvar[2]=1 Tvar[4]=3 Tvar[8]=8
        *       if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and          *       if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and 
        *  Tage[++cptcovage]=k         *  Tage[++cptcovage]=k
        *       if products, new covar are created after ncovcol with k1         *       if products, new covar are created after ncovcol + nqv (quanti fixed) with k1
        *  Tvar[k]=ncovcol+k1; # of the kth covariate product:  Tvar[5]=ncovcol+1=10  Tvar[6]=ncovcol+1=11         *  Tvar[k]=ncovcol+k1; # of the kth covariate product:  Tvar[5]=ncovcol+1=10  Tvar[6]=ncovcol+1=11
        *  Tprod[k1]=k; Tprod[1]=5 Tprod[2]= 6; gives the position of the k1th product         *  Tprod[k1]=k; Tprod[1]=5 Tprod[2]= 6; gives the position of the k1th product
        *  Tvard[k1][1]=m Tvard[k1][2]=m; Tvard[1][1]=5 (V5) Tvard[1][2]=6 Tvard[2][1]=7 (V7) Tvard[2][2]=8         *  Tvard[k1][1]=m Tvard[k1][2]=m; Tvard[1][1]=5 (V5) Tvard[1][2]=6 Tvard[2][1]=7 (V7) Tvard[2][2]=8
        *  Tvar[cptcovn+k2]=Tvard[k1][1];Tvar[cptcovn+k2+1]=Tvard[k1][2];         *  Tvar[cptcovn+k2]=Tvard[k1][1];Tvar[cptcovn+k2+1]=Tvard[k1][2];
        *  Tvar[8+1]=5;Tvar[8+2]=6;Tvar[8+3]=7;Tvar[8+4]=8 inverted         *  Tvar[8+1]=5;Tvar[8+2]=6;Tvar[8+3]=7;Tvar[8+4]=8 inverted
        *  V1   V2   V3   V4  V5  V6  V7  V8  V9  V10  V11         *  V1   V2   V3   V4  V5  V6  V7  V8  V9  V10  V11
        *  <          ncovcol=8                >         *  <          ncovcol=8  8 fixed covariate. Additional starts at 9 (V5*V6) and 10(V7*V8)              >
        *       Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8    d1   d1   d2  d2         *       Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8    d1   d1   d2  d2
        *          k=  1    2      3       4     5       6      7        8    9   10   11  12         *          k=  1    2      3       4     5       6      7        8    9   10   11  12
        *     Tvar[k]= 2    1      3       3    10      11      8        8    5    6    7   8         *     Tvard[k]= 2    1      3       3    10      11      8        8    5    6    7   8
        * p Tvar[1]@12={2,   1,     3,      3,  11,     10,     8,       8,   7,   8,   5,  6}         * p Tvar[1]@12={2,   1,     3,      3,   9,     10,     8,       8}
        * p Tprod[1]@2={                         6, 5}         * p Tprod[1]@2={                         6, 5}
        *p Tvard[1][1]@4= {7, 8, 5, 6}         *p Tvard[1][1]@4= {7, 8, 5, 6}
        * covar[k][i]= V2   V1      ?      V3    V5*V6?   V7*V8?  ?       V8            * covar[k][i]= V2   V1      ?      V3    V5*V6?   V7*V8?  ?       V8   
Line 11050  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 11308  Dummy[k] 0=dummy (0 1), 1 quantitative (
 Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product \n\  Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product \n\
 Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\  Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
 Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);  Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
   for(k=-1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;}    for(k=-1;k<=NCOVMAX; k++){ Fixed[k]=0; Dummy[k]=0;}
   for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0, ncovvt=0;k<=cptcovt; k++){ /* or cptocvt */    for(k=1;k<=NCOVMAX; k++){TvarFind[k]=0; TvarVind[k]=0;}
     for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0, ncovvt=0;k<=cptcovt; k++){ /* or cptocvt loop on k from model */
     if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */      if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
       Fixed[k]= 0;        Fixed[k]= 0;
       Dummy[k]= 0;        Dummy[k]= 0;
Line 11139  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 11398  Dummy[k] 0=dummy (0 1), 1 quantitative (
       TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */        TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */
       TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */        TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */
       /* Tmodeliqind[k]=nqtveff;/\* Only simple time varying quantitative variable *\/ */        /* Tmodeliqind[k]=nqtveff;/\* Only simple time varying quantitative variable *\/ */
       printf("Quasi TmodelQind[%d]=%d,Tvar[TmodelQind[%d]]=V%d, ncovcol=%d, nqv=%d, ntv=%d,Tvar[k]- ncovcol-nqv-ntv=%d\n",nqtveff,k,nqtveff,Tvar[k], ncovcol, nqv, ntv, Tvar[k]- ncovcol-nqv-ntv);        /* printf("Quasi TmodelQind[%d]=%d,Tvar[TmodelQind[%d]]=V%d, ncovcol=%d, nqv=%d, ntv=%d,Tvar[k]- ncovcol-nqv-ntv=%d\n",nqtveff,k,nqtveff,Tvar[k], ncovcol, nqv, ntv, Tvar[k]- ncovcol-nqv-ntv); */
       printf("Quasi TmodelInvQind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv-ntv);        /* printf("Quasi TmodelInvQind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv-ntv); */
     }else if (Typevar[k] == 1) {  /* product with age */      }else if (Typevar[k] == 1) {  /* product with age */
       ncova++;        ncova++;
       TvarA[ncova]=Tvar[k];        TvarA[ncova]=Tvar[k];
Line 11322  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 11581  Dummy[k] 0=dummy (0 1), 1 quantitative (
       printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);        printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
       fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);        fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
     }      }
     printf("Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[k],Dummy[k]);      /* printf("Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[k],Dummy[k]); */
     printf("           modell[%d].maintype=%d, modell[%d].subtype=%d\n",k,modell[k].maintype,k,modell[k].subtype);      /* printf("           modell[%d].maintype=%d, modell[%d].subtype=%d\n",k,modell[k].maintype,k,modell[k].subtype); */
     fprintf(ficlog,"Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[k],Dummy[k]);      fprintf(ficlog,"Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[k],Dummy[k]);
   }    }
   /* Searching for doublons in the model */    /* Searching for doublons in the model */
Line 12477  int main(int argc, char *argv[]) Line 12736  int main(int argc, char *argv[])
     numlinepar++;      numlinepar++;
     if(line[1]=='q'){ /* This #q will quit imach (the answer is q) */      if(line[1]=='q'){ /* This #q will quit imach (the answer is q) */
       z[0]=line[1];        z[0]=line[1];
       }else if(line[1]=='d'){ /* For debugging individual values of covariates in ficresilk */
         debugILK=1;printf("DebugILK\n");
     }      }
     /* printf("****line [1] = %c \n",line[1]); */      /* printf("****line [1] = %c \n",line[1]); */
     fputs(line, stdout);      fputs(line, stdout);
Line 12490  int main(int argc, char *argv[]) Line 12751  int main(int argc, char *argv[])
   covar=matrix(0,NCOVMAX,firstobs,lastobs);  /**< used in readdata */    covar=matrix(0,NCOVMAX,firstobs,lastobs);  /**< used in readdata */
   if(nqv>=1)coqvar=matrix(1,nqv,firstobs,lastobs);  /**< Fixed quantitative covariate */    if(nqv>=1)coqvar=matrix(1,nqv,firstobs,lastobs);  /**< Fixed quantitative covariate */
   if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,firstobs,lastobs);  /**< Time varying quantitative covariate */    if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,firstobs,lastobs);  /**< Time varying quantitative covariate */
   if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,firstobs,lastobs);  /**< Time varying covariate (dummy and quantitative)*/    /* if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,firstobs,lastobs);  /\**< Time varying covariate (dummy and quantitative)*\/ */
     if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,ncovcol+nqv+1,ncovcol+nqv+ntv+nqtv,firstobs,lastobs);  /**< Might be better */
   cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/    cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/
   /* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5    /* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5
      v1+v2*age+v2*v3 makes cptcovn = 3       v1+v2*age+v2*v3 makes cptcovn = 3
Line 13774  Please run with mle=-1 to get a correct Line 14036  Please run with mle=-1 to get a correct
       case 13:        case 13:
         num_filled=sscanf(line,"result:%[^\n]\n",resultlineori);          num_filled=sscanf(line,"result:%[^\n]\n",resultlineori);
         nresult++; /* Sum of resultlines */          nresult++; /* Sum of resultlines */
         printf("Result %d: result:%s\n",nresult, resultlineori);          /* printf("Result %d: result:%s\n",nresult, resultlineori); */
         /* removefirstspace(&resultlineori); */          /* removefirstspace(&resultlineori); */
                   
         if(strstr(resultlineori,"v") !=0){          if(strstr(resultlineori,"v") !=0){
Line 13783  Please run with mle=-1 to get a correct Line 14045  Please run with mle=-1 to get a correct
           return 1;            return 1;
         }          }
         trimbb(resultline, resultlineori); /* Suppressing double blank in the resultline */          trimbb(resultline, resultlineori); /* Suppressing double blank in the resultline */
         printf("Decoderesult resultline=\"%s\" resultlineori=\"%s\"\n", resultline, resultlineori);          /* printf("Decoderesult resultline=\"%s\" resultlineori=\"%s\"\n", resultline, resultlineori); */
         if(nresult > MAXRESULTLINESPONE-1){          if(nresult > MAXRESULTLINESPONE-1){
           printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINESPONE-1,nresult,rfileres);            printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINESPONE-1,nresult,rfileres);
           fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINESPONE-1,nresult,rfileres);            fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINESPONE-1,nresult,rfileres);
Line 14120  Please run with mle=-1 to get a correct Line 14382  Please run with mle=-1 to get a correct
         /* Tvresult[nres][j] Name of the variable at position j in this resultline */          /* Tvresult[nres][j] Name of the variable at position j in this resultline */
         /* Tresult[nres][j] Value of this variable at position j could be a float if quantitative  */          /* Tresult[nres][j] Value of this variable at position j could be a float if quantitative  */
 /* We give up with the combinations!! */  /* We give up with the combinations!! */
         printf("\n j=%d In computing T_ Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=%d cptcovs=%d, cptcoveff=%d Fixed[modelresult[nres][j]]=%d\n", j, nres, j, Dummy[modelresult[nres][j]],nres,j,modelresult[nres][j],cptcovs, cptcoveff,Fixed[modelresult[nres][j]]);  /* end if dummy  or quanti */          /* if(debugILK) */
           /*   printf("\n j=%d In computing T_ Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=%d cptcovs=%d, cptcoveff=%d Fixed[modelresult[nres][j]]=%d\n", j, nres, j, Dummy[modelresult[nres][j]],nres,j,modelresult[nres][j],cptcovs, cptcoveff,Fixed[modelresult[nres][j]]);  /\* end if dummy  or quanti *\/ */
   
         if(Dummy[modelresult[nres][j]]==0){/* Dummy variable of the variable in position modelresult in the model corresponding to j in resultline  */          if(Dummy[modelresult[nres][j]]==0){/* Dummy variable of the variable in position modelresult in the model corresponding to j in resultline  */
           printf("V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */            /* printf("V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][j]); /\* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  *\/ */ /* TinvDoQresult[nres][Name of the variable] */
           fprintf(ficlog,"V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */            printf("V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); /* Output of each value for the combination TKresult[nres], ordered by the covariate values in the resultline  */
           fprintf(ficrest,"V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */            fprintf(ficlog,"V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
             fprintf(ficrest,"V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
           if(Fixed[modelresult[nres][j]]==0){ /* Fixed */            if(Fixed[modelresult[nres][j]]==0){ /* Fixed */
             printf("fixed ");fprintf(ficlog,"fixed ");fprintf(ficrest,"fixed ");              printf("fixed ");fprintf(ficlog,"fixed ");fprintf(ficrest,"fixed ");
           }else{            }else{
Line 14294  Please run with mle=-1 to get a correct Line 14558  Please run with mle=-1 to get a correct
   free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);    free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);    free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);    free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
   if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,firstobs,lastobs);    /* if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,firstobs,lastobs); */
     if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,ncovcol+nqv+1,ncovcol+nqv+ntv+nqtv,firstobs,lastobs);
   if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,firstobs,lastobs);    if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,firstobs,lastobs);
   if(nqv>=1)free_matrix(coqvar,1,nqv,firstobs,lastobs);    if(nqv>=1)free_matrix(coqvar,1,nqv,firstobs,lastobs);
   free_matrix(covar,0,NCOVMAX,firstobs,lastobs);    free_matrix(covar,0,NCOVMAX,firstobs,lastobs);

Removed from v.1.339  
changed lines
  Added in v.1.345


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