Diff for /imach/src/imach.c between versions 1.333 and 1.336

version 1.333, 2022/08/21 09:10:30 version 1.336, 2022/08/31 09:52:36
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.336  2022/08/31 09:52:36  brouard
     *** empty log message ***
   
     Revision 1.335  2022/08/31 08:23:16  brouard
     Summary: improvements...
   
     Revision 1.334  2022/08/25 09:08:41  brouard
     Summary: In progress for quantitative
   
   Revision 1.333  2022/08/21 09:10:30  brouard    Revision 1.333  2022/08/21 09:10:30  brouard
   * src/imach.c (Module): Version 0.99r33 A lot of changes in    * src/imach.c (Module): Version 0.99r33 A lot of changes in
   reassigning covariates: my first idea was that people will always    reassigning covariates: my first idea was that people will always
Line 1290  int erreur=0, nberr=0, nbwarn=0; /* Erro Line 1299  int erreur=0, nberr=0, nbwarn=0; /* Erro
 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 */
 int cptcovn=0; /**< cptcovn decodemodel: number of covariates k of the models excluding age*products =6 and age*age */  int cptcovn=0; /**< cptcovn decodemodel: number of covariates k of the models excluding age*products =6 and age*age but including products */
 int cptcovt=0; /**< cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */  int cptcovt=0; /**< cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */
 int cptcovs=0; /**< cptcovs number of simple covariates in the model V2+V1 =2 (dummy or quantit or time varying) */  int cptcovs=0; /**< cptcovs number of SIMPLE covariates in the model V2+V1 =2 (dummy or quantit or time varying) */
 int cptcovsnq=0; /**< cptcovsnq number of simple covariates in the model but non quantitative V2+V1 =2 */  int cptcovsnq=0; /**< cptcovsnq number of SIMPLE covariates in the model but non quantitative V2+V1 =2 */
 int cptcovage=0; /**< Number of covariates with age: V3*age only =1 */  int cptcovage=0; /**< Number of covariates with age: V3*age only =1 */
 int cptcovprodnoage=0; /**< Number of covariate products without age */     int cptcovprodnoage=0; /**< Number of covariate products without age */   
 int cptcoveff=0; /* Total number of covariates to vary for printing results (2**cptcoveff combinations of dummies)(computed in tricode as cptcov) */  int cptcoveff=0; /* Total number of single dummy covariates (fixed or time varying) to vary for printing results (2**cptcoveff combinations of dummies)(computed in tricode as cptcov) */
 int ncovf=0; /* Total number of effective fixed covariates (dummy or quantitative) in the model */  int ncovf=0; /* Total number of effective fixed covariates (dummy or quantitative) in the model */
 int ncovv=0; /* Total number of effective (wave) varying covariates (dummy or quantitative) in the model */  int ncovv=0; /* Total number of effective (wave) varying covariates (dummy or quantitative) in the model */
 int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (dummy of quantitative) in the model */  int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (dummy of quantitative) in the model */
Line 1307  int nqfveff=0; /**< nqfveff Number of Qu Line 1316  int nqfveff=0; /**< nqfveff Number of Qu
 int ntveff=0; /**< ntveff number of effective time varying variables */  int ntveff=0; /**< ntveff number of effective time varying variables */
 int nqtveff=0; /**< ntqveff number of effective time varying quantitative variables */  int nqtveff=0; /**< ntqveff number of effective time varying quantitative variables */
 int cptcov=0; /* Working variable */  int cptcov=0; /* Working variable */
   int firstobs=1, lastobs=10; /* nobs = lastobs-firstobs+1 declared globally ;*/
 int nobs=10;  /* Number of observations in the data lastobs-firstobs */  int nobs=10;  /* Number of observations in the data lastobs-firstobs */
 int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */  int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */
 int npar=NPARMAX; /* Number of parameters (nlstate+ndeath-1)*nlstate*ncovmodel; */  int npar=NPARMAX; /* Number of parameters (nlstate+ndeath-1)*nlstate*ncovmodel; */
Line 1524  int *TvarsQind; Line 1534  int *TvarsQind;
 #define MAXRESULTLINESPONE 10+1  #define MAXRESULTLINESPONE 10+1
 int nresult=0;  int nresult=0;
 int parameterline=0; /* # of the parameter (type) line */  int parameterline=0; /* # of the parameter (type) line */
 int TKresult[MAXRESULTLINESPONE];  int TKresult[MAXRESULTLINESPONE]; /* TKresult[nres]=k for each resultline nres give the corresponding combination of dummies */
 int resultmodel[MAXRESULTLINESPONE][NCOVMAX];/* resultmodel[k1]=k3: k1th position in the model correspond to the k3 position in the resultline */  int resultmodel[MAXRESULTLINESPONE][NCOVMAX];/* resultmodel[k1]=k3: k1th position in the model corresponds to the k3 position in the resultline */
 int Tresult[MAXRESULTLINESPONE][NCOVMAX];/* For dummy variable , value (output) */  int modelresult[MAXRESULTLINESPONE][NCOVMAX];/* modelresult[k3]=k1: k1th position in the model corresponds to the k3 position in the resultline */
   int Tresult[MAXRESULTLINESPONE][NCOVMAX];/* Tresult[nres][result_position]= value of the dummy variable at the result_position in the nres resultline */
 int Tinvresult[MAXRESULTLINESPONE][NCOVMAX];/* Tinvresult[nres][Name of a dummy variable]= value of the variable in the result line  */  int Tinvresult[MAXRESULTLINESPONE][NCOVMAX];/* Tinvresult[nres][Name of a dummy variable]= value of the variable in the result line  */
 double TinvDoQresult[MAXRESULTLINESPONE][NCOVMAX];/* TinvDoQresult[nres][Name of a Dummy or Q variable]= value of the variable in the result line */  double TinvDoQresult[MAXRESULTLINESPONE][NCOVMAX];/* TinvDoQresult[nres][Name of a Dummy or Q variable]= value of the variable in the result line */
 int Tvresult[MAXRESULTLINESPONE][NCOVMAX]; /* For dummy variable , variable # (output) */  int Tvresult[MAXRESULTLINESPONE][NCOVMAX]; /* Tvresult[nres][result_position]= name of the dummy variable at the result_position in the nres resultline */
 double Tqresult[MAXRESULTLINESPONE][NCOVMAX]; /* Tqresult[nres][result_position]= value of the variable at the result_position in the nres resultline */  double Tqresult[MAXRESULTLINESPONE][NCOVMAX]; /* Tqresult[nres][result_position]= value of the variable at the result_position in the nres resultline */
 double Tqinvresult[MAXRESULTLINESPONE][NCOVMAX]; /* For quantitative variable , value (output) */  double Tqinvresult[MAXRESULTLINESPONE][NCOVMAX]; /* For quantitative variable , value (output) */
 int Tvqresult[MAXRESULTLINESPONE][NCOVMAX]; /* Tvqresult[nres][result_position]= id of the variable at the result_position in the nres resultline */  int Tvqresult[MAXRESULTLINESPONE][NCOVMAX]; /* Tvqresult[nres][result_position]= id of the variable at the result_position in the nres resultline */
Line 3500  double **matprod2(double **out, double * Line 3511  double **matprod2(double **out, double *
   
 double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij, int nres )  double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij, int nres )
 {  {
   /* Computes the transition matrix starting at age 'age' and dummies values in each resultline (loop on ij to find the corresponding combination) to over     /* Already optimized with precov.
        Computes the transition matrix starting at age 'age' and dummies values in each resultline (loop on ij to find the corresponding combination) to over 
      'nhstepm*hstepm*stepm' months (i.e. until       'nhstepm*hstepm*stepm' months (i.e. until
      age (in years)  age+nhstepm*hstepm*stepm/12) by multiplying        age (in years)  age+nhstepm*hstepm*stepm/12) by multiplying 
      nhstepm*hstepm matrices.        nhstepm*hstepm matrices. 
Line 3831  double ***hbxij(double ***po, int nhstep Line 3843  double ***hbxij(double ***po, int nhstep
 /*************** log-likelihood *************/  /*************** log-likelihood *************/
 double func( double *x)  double func( double *x)
 {  {
   int i, ii, j, k, mi, d, kk;    int i, ii, j, k, mi, d, kk, kf=0;
   int ioffset=0;    int ioffset=0;
   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 */
   int s1, s2;    int s1, s2;
   int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate, quantitative time varying covariate */    int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate, quantitative time varying covariate */
   
   double bbh, survp;    double bbh, survp;
   long ipmx;  
   double agexact;    double agexact;
     double agebegin, ageend;
   /*extern weight */    /*extern weight */
   /* We are differentiating ll according to initial status */    /* We are differentiating ll according to initial status */
   /*  for (i=1;i<=npar;i++) printf("%f ", x[i]);*/    /*  for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
Line 3863  double func( double *x) Line 3876  double func( double *x)
       */        */
       ioffset=2+nagesqr ;        ioffset=2+nagesqr ;
    /* Fixed */     /* Fixed */
       for (k=1; k<=ncovf;k++){ /* For each fixed covariate dummu or quant or prod */        for (kf=1; kf<=ncovf;kf++){ /* For each fixed covariate dummu 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 */
         /* TvarFind;  TvarFind[1]=6,  TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod)  */          /* TvarFind;  TvarFind[1]=6,  TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod)  */
         cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (TvarFind[1]=6)*/          cov[ioffset+TvarFind[kf]]=covar[Tvar[TvarFind[kf]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (TvarFind[1]=6)*/
         /* V1*V2 (7)  TvarFind[2]=7, TvarFind[3]=9 */          /* V1*V2 (7)  TvarFind[2]=7, TvarFind[3]=9 */
       }        }
       /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4]         /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] 
Line 3883  double func( double *x) Line 3896  double func( double *x)
          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]
       */        */
       for(mi=1; mi<= wav[i]-1; mi++){        for(mi=1; mi<= wav[i]-1; mi++){  /* Varying with waves */
         /* Wave varying (but not age varying) */
         for(k=1; k <= ncovv ; k++){ /* Varying  covariates in the model (single and product but no age )"V5+V4+V3+V4*V3+V5*age+V1*age+V1" +TvarVind 1,2,3,4(V4*V3)  Tvar[1]@7{5, 4, 3, 6, 5, 1, 1 ; 6 because the created covar is after V5 and is 6, minus 1+1, 3,2,1,4 positions in cotvar*/          for(k=1; k <= ncovv ; k++){ /* Varying  covariates in the model (single and product but no age )"V5+V4+V3+V4*V3+V5*age+V1*age+V1" +TvarVind 1,2,3,4(V4*V3)  Tvar[1]@7{5, 4, 3, 6, 5, 1, 1 ; 6 because the created covar is after V5 and is 6, minus 1+1, 3,2,1,4 positions in cotvar*/
           /* 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];
Line 3893  double func( double *x) Line 3907  double func( double *x)
             oldm[ii][j]=(ii==j ? 1.0 : 0.0);              oldm[ii][j]=(ii==j ? 1.0 : 0.0);
             savm[ii][j]=(ii==j ? 1.0 : 0.0);              savm[ii][j]=(ii==j ? 1.0 : 0.0);
           }            }
   
           agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */
           ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */
         for(d=0; d<dh[mi][i]; d++){          for(d=0; d<dh[mi][i]; d++){
           newm=savm;            newm=savm;
           agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;            agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
Line 3983  double func( double *x) Line 4000  double func( double *x)
           /*survp += out[s1][j]; */            /*survp += out[s1][j]; */
           lli= log(survp);            lli= log(survp);
         }          }
         else if  (s2==-4) {           /* else if  (s2==-4) {  */
           for (j=3,survp=0. ; j<=nlstate; j++)            /*   for (j=3,survp=0. ; j<=nlstate; j++)   */
             survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];          /*     survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; */
           lli= log(survp);           /*   lli= log(survp);  */
         }           /* }  */
         else if  (s2==-5) {           /* else if  (s2==-5) {  */
           for (j=1,survp=0. ; j<=2; j++)            /*   for (j=1,survp=0. ; j<=2; j++)   */
             survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];          /*     survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; */
           lli= log(survp);           /*   lli= log(survp);  */
         }           /* }  */
         else{          else{
           lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */            lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
           /*  lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */            /*  lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */
Line 4167  double func( double *x) Line 4184  double func( double *x)
 double funcone( double *x)  double funcone( double *x)
 {  {
   /* Same as func but slower because of a lot of printf and if */    /* Same as func but slower because of a lot of printf and if */
   int i, ii, j, k, mi, d, kk;    int i, ii, j, k, mi, d, kk, kf=0;
   int ioffset=0;    int ioffset=0;
   double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];    double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
   double **out;    double **out;
Line 4190  double funcone( double *x) Line 4207  double funcone( double *x)
   for(k=1; k<=nlstate; k++) ll[k]=0.;    for(k=1; k<=nlstate; k++) ll[k]=0.;
   ioffset=0;    ioffset=0;
   for (i=1,ipmx=0, sw=0.; i<=imx; i++){    for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       /* Computes the values of the ncovmodel covariates of the model
          depending if the covariates are fixed or varying (age dependent) and stores them in cov[]
          Then computes with function pmij which return a matrix p[i][j] giving the elementary probability
          to be observed in j being in i according to the model.
       */
     /* ioffset=2+nagesqr+cptcovage; */      /* ioffset=2+nagesqr+cptcovage; */
     ioffset=2+nagesqr;      ioffset=2+nagesqr;
     /* Fixed */      /* Fixed */
     /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */      /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */
     /* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */      /* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */
     for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */      for (kf=1; kf<=ncovf;kf++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */
       cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/        cov[ioffset+TvarFind[kf]]=covar[Tvar[TvarFind[kf]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/
 /*    cov[ioffset+TvarFind[1]]=covar[Tvar[TvarFind[1]]][i];  */  /*    cov[ioffset+TvarFind[1]]=covar[Tvar[TvarFind[1]]][i];  */
 /*    cov[2+6]=covar[Tvar[6]][i];  */  /*    cov[2+6]=covar[Tvar[6]][i];  */
 /*    cov[2+6]=covar[2][i]; V2  */  /*    cov[2+6]=covar[2][i]; V2  */
Line 4207  double funcone( double *x) Line 4229  double funcone( double *x)
 /*    cov[2+9]=covar[Tvar[9]][i];  */  /*    cov[2+9]=covar[Tvar[9]][i];  */
 /*    cov[2+9]=covar[1][i]; V1  */  /*    cov[2+9]=covar[1][i]; V1  */
     }      }
         /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] 
            is 5, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]=6 
            has been calculated etc */
         /* For an individual i, wav[i] gives the number of effective waves */
         /* We compute the contribution to Likelihood of each effective transition
            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];
            s2=s[mw[mi+1][i]][i];
            And the iv th varying covariate is the cotvar[mw[mi+1][i]][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 */
     /* 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 *\/ */
     /*   cov[++ioffset]=coqvar[TvarFQ[k]][i];/\* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V2 and V1*V2 is fixed (k=6 and 7?)*\/ */      /*   cov[++ioffset]=coqvar[TvarFQ[k]][i];/\* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V2 and V1*V2 is fixed (k=6 and 7?)*\/ */
     /* } */      /* } */
Line 4264  double funcone( double *x) Line 4299  double funcone( double *x)
         savm=oldm;          savm=oldm;
         oldm=newm;          oldm=newm;
       } /* end mult */        } /* end mult */
                 /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
           /* But now since version 0.9 we anticipate for bias at large stepm.
            * If stepm is larger than one month (smallest stepm) and if the exact delay 
            * (in months) between two waves is not a multiple of stepm, we rounded to 
            * the nearest (and in case of equal distance, to the lowest) interval but now
            * we keep into memory the bias bh[mi][i] and also the previous matrix product
            * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the
            * probability in order to take into account the bias as a fraction of the way
                                    * from savm to out if bh is negative or even beyond if bh is positive. bh varies
                                    * -stepm/2 to stepm/2 .
                                    * For stepm=1 the results are the same as for previous versions of Imach.
                                    * For stepm > 1 the results are less biased than in previous versions. 
                                    */
       s1=s[mw[mi][i]][i];        s1=s[mw[mi][i]][i];
       s2=s[mw[mi+1][i]][i];        s2=s[mw[mi+1][i]][i];
       /* if(s2==-1){ */        /* if(s2==-1){ */
Line 4296  double funcone( double *x) Line 4343  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("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("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])); */
       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\ */
    /* %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, */
    /*             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); */
         }          }
         fprintf(ficresilk," %10.6f\n", -llt);          fprintf(ficresilk," %10.6f\n", -llt);
           /* printf(" %10.6f\n", -llt); */
       }        }
         } /* 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];
 /* printf("l1=%f l2=%f ",ll[1],ll[2]); */  /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
 l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */    l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
 if(globpr==0){ /* First time we count the contributions and weights */    if(globpr==0){ /* First time we count the contributions and weights */
         gipmx=ipmx;      gipmx=ipmx;
         gsw=sw;      gsw=sw;
 }    }
 return -l;  return -l;
 }  }
   
Line 4975  Title=%s <br>Datafile=%s Firstpass=%d La Line 5028  Title=%s <br>Datafile=%s Firstpass=%d La
   j1=0;    j1=0;
       
   /* j=ncoveff;  /\* Only fixed dummy covariates *\/ */    /* j=ncoveff;  /\* Only fixed dummy covariates *\/ */
   j=cptcoveff;  /* Only dummy covariates of the model */    j=cptcoveff;  /* Only simple dummy covariates used in the model */
   /* j=cptcovn;  /\* Only dummy covariates of the model *\/ */    /* j=cptcovn;  /\* Only dummy covariates of the model *\/ */
   if (cptcovn<1) {j=1;ncodemax[1]=1;}    if (cptcovn<1) {j=1;ncodemax[1]=1;}
       
Line 4996  Title=%s <br>Datafile=%s Firstpass=%d La Line 5049  Title=%s <br>Datafile=%s Firstpass=%d La
   
   /* if a constant only model, one pass to compute frequency tables and to write it on ficresp */    /* if a constant only model, one pass to compute frequency tables and to write it on ficresp */
   /* Loop on nj=1 or 2 if dummy covariates j!=0    /* Loop on nj=1 or 2 if dummy covariates j!=0
    *   Loop on j1(1 to 2**cptcovn) covariate combination     *   Loop on j1(1 to 2**cptcoveff) covariate combination
    *     freq[s1][s2][iage] =0.     *     freq[s1][s2][iage] =0.
    *     Loop on iind     *     Loop on iind
    *       ++freq[s1][s2][iage] weighted     *       ++freq[s1][s2][iage] weighted
Line 5021  Title=%s <br>Datafile=%s Firstpass=%d La Line 5074  Title=%s <br>Datafile=%s Firstpass=%d La
     if(nj==1)      if(nj==1)
       j=0;  /* First pass for the constant */        j=0;  /* First pass for the constant */
     else{      else{
       j=cptcovs; /* Other passes for the covariate values */        j=cptcoveff; /* Other passes for the covariate values number of simple covariates in the model V2+V1 =2 (simple dummy fixed or time varying) */
     }      }
     first=1;      first=1;
     for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on all dummy covariates combination of the model, ie excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */      for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on all dummy covariates combination of the model, ie excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */
Line 5057  Title=%s <br>Datafile=%s Firstpass=%d La Line 5110  Title=%s <br>Datafile=%s Firstpass=%d La
         bool=1;          bool=1;
         if(j !=0){          if(j !=0){
           if(anyvaryingduminmodel==0){ /* If All fixed covariates */            if(anyvaryingduminmodel==0){ /* If All fixed covariates */
             if (cptcovn >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */              if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
               for (z1=1; z1<=cptcovn; z1++) { /* loops on covariates in the model */                for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */
                 /* if(Tvaraff[z1] ==-20){ */                  /* if(Tvaraff[z1] ==-20){ */
                 /*       /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */                  /*       /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
                 /* }else  if(Tvaraff[z1] ==-10){ */                  /* }else  if(Tvaraff[z1] ==-10){ */
                 /*       /\* sumnew+=coqvar[z1][iind]; *\/ */                  /*       /\* sumnew+=coqvar[z1][iind]; *\/ */
                 /* }else  */ /* TODO TODO codtabm(j1,z1) or codtabm(j1,Tvaraff[z1]]z1)*/                  /* }else  */ /* TODO TODO codtabm(j1,z1) or codtabm(j1,Tvaraff[z1]]z1)*/
                   /* if( iind >=imx-3) printf("Searching error iind=%d Tvaraff[z1]=%d covar[Tvaraff[z1]][iind]=%.f TnsdVar[Tvaraff[z1]]=%d, cptcoveff=%d, cptcovs=%d \n",iind, Tvaraff[z1], covar[Tvaraff[z1]][iind],TnsdVar[Tvaraff[z1]],cptcoveff, cptcovs); */
                   if(Tvaraff[z1]<1 || Tvaraff[z1]>=NCOVMAX)
                     printf("Error Tvaraff[z1]=%d<1 or >=%d, cptcoveff=%d model=%s\n",Tvaraff[z1],NCOVMAX, cptcoveff, model);
                 if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]){ /* for combination j1 of covariates */                  if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]){ /* for combination j1 of covariates */
                   /* Tests if the value of the covariate z1 for this individual iind responded to combination j1 (V4=1 V3=0) */                    /* Tests if the value of the covariate z1 for this individual iind responded to combination j1 (V4=1 V3=0) */
                   bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */                    bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */
Line 5073  Title=%s <br>Datafile=%s Firstpass=%d La Line 5129  Title=%s <br>Datafile=%s Firstpass=%d La
                   /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/                    /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
                 } /* Onlyf fixed */                  } /* Onlyf fixed */
               } /* end z1 */                } /* end z1 */
             } /* cptcovn > 0 */              } /* cptcoveff > 0 */
           } /* end any */            } /* end any */
         }/* end j==0 */          }/* end j==0 */
         if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */          if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */
Line 5082  Title=%s <br>Datafile=%s Firstpass=%d La Line 5138  Title=%s <br>Datafile=%s Firstpass=%d La
             m=mw[mi][iind];              m=mw[mi][iind];
             if(j!=0){              if(j!=0){
               if(anyvaryingduminmodel==1){ /* Some are varying covariates */                if(anyvaryingduminmodel==1){ /* Some are varying covariates */
                 for (z1=1; z1<=cptcovn; 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;
                     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 
Line 5090  Title=%s <br>Datafile=%s Firstpass=%d La Line 5146  Title=%s <br>Datafile=%s Firstpass=%d La
                                                                                       constant and age model which counts them. */                                                                                        constant and age model which counts them. */
                       bool=0; /* not selected */                        bool=0; /* not selected */
                   }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */                    }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */
                     if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) {                      /* i1=Tvaraff[z1]; */
                       /* i2=TnsdVar[i1]; */
                       /* i3=nbcode[i1][i2]; */
                       /* i4=covar[i1][iind]; */
                       /* if(i4 != i3){ */
                       if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) { /* Bug valgrind */
                       bool=0;                        bool=0;
                     }                      }
                   }                    }
Line 5151  Title=%s <br>Datafile=%s Firstpass=%d La Line 5212  Title=%s <br>Datafile=%s Firstpass=%d La
               
               
       /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/        /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
       if(cptcovn==0 && nj==1) /* no covariate and first pass */        if(cptcoveff==0 && nj==1) /* no covariate and first pass */
         pstamp(ficresp);          pstamp(ficresp);
       if  (cptcovn>0 && j!=0){        if  (cptcoveff>0 && j!=0){
         pstamp(ficresp);          pstamp(ficresp);
         printf( "\n#********** Variable ");           printf( "\n#********** Variable "); 
         fprintf(ficresp, "\n#********** Variable ");           fprintf(ficresp, "\n#********** Variable "); 
Line 5206  Title=%s <br>Datafile=%s Firstpass=%d La Line 5267  Title=%s <br>Datafile=%s Firstpass=%d La
       /* } */        /* } */
   
       fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");        fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
       if((cptcovn==0 && nj==1)|| nj==2 ) /* no covariate and first pass */        if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */
         fprintf(ficresp, " Age");          fprintf(ficresp, " Age");
       if(nj==2) for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " V%d=%d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);        if(nj==2) for (z1=1; z1<=cptcoveff; z1++) {
             printf(" V%d=%d, z1=%d, Tvaraff[z1]=%d, j1=%d, TnsdVar[Tvaraff[%d]]=%d |",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])], z1, Tvaraff[z1], j1,z1,TnsdVar[Tvaraff[z1]]);
             fprintf(ficresp, " V%d=%d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
           }
       for(i=1; i<=nlstate;i++) {        for(i=1; i<=nlstate;i++) {
         if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp," Prev(%d)  N(%d)  N  ",i,i);          if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," Prev(%d)  N(%d)  N  ",i,i);
         fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);          fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);
       }        }
       if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n");        if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n");
       fprintf(ficresphtm, "\n");        fprintf(ficresphtm, "\n");
               
       /* Header of frequency table by age */        /* Header of frequency table by age */
Line 5281  Title=%s <br>Datafile=%s Firstpass=%d La Line 5345  Title=%s <br>Datafile=%s Firstpass=%d La
         }          }
                   
         /* Writing ficresp */          /* Writing ficresp */
         if(cptcovn==0 && nj==1){ /* no covariate and first pass */          if(cptcoveff==0 && nj==1){ /* no covariate and first pass */
           if( iage <= iagemax){            if( iage <= iagemax){
             fprintf(ficresp," %d",iage);              fprintf(ficresp," %d",iage);
           }            }
         }else if( nj==2){          }else if( nj==2){
           if( iage <= iagemax){            if( iage <= iagemax){
             fprintf(ficresp," %d",iage);              fprintf(ficresp," %d",iage);
             for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);              for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
           }            }
         }          }
         for(s1=1; s1 <=nlstate ; s1++){          for(s1=1; s1 <=nlstate ; s1++){
Line 5303  Title=%s <br>Datafile=%s Firstpass=%d La Line 5367  Title=%s <br>Datafile=%s Firstpass=%d La
           }            }
           if( iage <= iagemax){            if( iage <= iagemax){
             if(pos>=1.e-5){              if(pos>=1.e-5){
               if(cptcovn==0 && nj==1){ /* no covariate and first pass */                if(cptcoveff==0 && nj==1){ /* no covariate and first pass */
                 fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta);                  fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta);
               }else if( nj==2){                }else if( nj==2){
                 fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta);                  fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta);
Line 5312  Title=%s <br>Datafile=%s Firstpass=%d La Line 5376  Title=%s <br>Datafile=%s Firstpass=%d La
               /*probs[iage][s1][j1]= pp[s1]/pos;*/                /*probs[iage][s1][j1]= pp[s1]/pos;*/
               /*printf("\niage=%d s1=%d j1=%d %.5f %.0f %.0f %f",iage,s1,j1,pp[s1]/pos, pp[s1],pos,probs[iage][s1][j1]);*/                /*printf("\niage=%d s1=%d j1=%d %.5f %.0f %.0f %f",iage,s1,j1,pp[s1]/pos, pp[s1],pos,probs[iage][s1][j1]);*/
             } else{              } else{
               if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta);                if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta);
               fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[s1][iage],pospropta);                fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[s1][iage],pospropta);
             }              }
           }            }
Line 5338  Title=%s <br>Datafile=%s Firstpass=%d La Line 5402  Title=%s <br>Datafile=%s Firstpass=%d La
         }          }
         fprintf(ficresphtmfr,"</tr>\n ");          fprintf(ficresphtmfr,"</tr>\n ");
         fprintf(ficresphtm,"</tr>\n");          fprintf(ficresphtm,"</tr>\n");
         if((cptcovn==0 && nj==1)|| nj==2 ) {          if((cptcoveff==0 && nj==1)|| nj==2 ) {
           if(iage <= iagemax)            if(iage <= iagemax)
             fprintf(ficresp,"\n");              fprintf(ficresp,"\n");
         }          }
Line 5596  void prevalence(double ***probs, double Line 5660  void prevalence(double ***probs, double
   if (cptcovn<1) {j=1;ncodemax[1]=1;}    if (cptcovn<1) {j=1;ncodemax[1]=1;}
       
   first=0;    first=0;
   for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */    for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of simple dummy covariates */
     for (i=1; i<=nlstate; i++)        for (i=1; i<=nlstate; i++)  
       for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++)        for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++)
         prop[i][iage]=0.0;          prop[i][iage]=0.0;
Line 5920  void  concatwav(int wav[], int **dh, int Line 5984  void  concatwav(int wav[], int **dh, int
        nbcode[k][j]=0; /* Valgrind */         nbcode[k][j]=0; /* Valgrind */
   
    /* 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++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */     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;
      if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */        if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ 
        switch(Fixed[k]) {         switch(Fixed[k]) {
Line 6018  void  concatwav(int wav[], int **dh, int Line 6082  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){ /* Dummy covariate and not age product */        if(Dummy[k]==1 && Typevar[k] !=1){ /* 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){
              printf("Error k=%d \n",k);
              exit(1);
            }
          if(isnan(covar[Tvar[k]][i])){           if(isnan(covar[Tvar[k]][i])){
            printf("ERROR, IMaCh doesn't treat fixed quantitative covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i);             printf("ERROR, IMaCh doesn't treat fixed quantitative covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i);
            fprintf(ficlog,"ERROR, currently IMaCh doesn't treat covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i);             fprintf(ficlog,"ERROR, currently IMaCh doesn't treat covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i);
Line 6027  void  concatwav(int wav[], int **dh, int Line 6095  void  concatwav(int wav[], int **dh, int
            exit(1);             exit(1);
          }           }
        }         }
      }       } /* end Quanti */
    } /* end of loop on model-covariate k. nbcode[Tvark][1]=-1, nbcode[Tvark][1]=0 and nbcode[Tvark][2]=1 sets the value of covariate k*/       } /* end of loop on model-covariate k. nbcode[Tvark][1]=-1, nbcode[Tvark][1]=0 and nbcode[Tvark][2]=1 sets the value of covariate k*/  
       
    for (k=-1; k< maxncov; k++) Ndum[k]=0;      for (k=-1; k< maxncov; k++) Ndum[k]=0; 
Line 6041  void  concatwav(int wav[], int **dh, int Line 6109  void  concatwav(int wav[], int **dh, int
       
    ij=0;     ij=0;
    /* for (i=0; i<=  maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) *\/ */     /* for (i=0; i<=  maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) *\/ */
    for (k=1; k<=  cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */     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 */
        /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
      /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/       /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
      /* if((Ndum[i]!=0) && (i<=ncovcol)){  /\* Tvar[i] <= ncovmodel ? *\/ */       /* if((Ndum[i]!=0) && (i<=ncovcol)){  /\* Tvar[i] <= ncovmodel ? *\/ */
      if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){  /* Only Dummy and non empty in the model */       if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){  /* Only Dummy simple and non empty in the model */
          /* Typevar[k] =0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */
          /* Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product*/
        /* If product not in single variable we don't print results */         /* If product not in single variable we don't print results */
        /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/         /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
        ++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */         ++ij;/*    V5 + V4 + V3 + V4*V3 + V5*age + V2 +  V1*V2 + V1*age + V1, *//* V5 quanti, V2 quanti, V4, V3, V1 dummies */
          /* k=       1    2   3     4       5       6      7       8        9  */
          /* Tvar[k]= 5    4    3    6       5       2      7       1        1  */
          /* ij            1    2                                            3  */  
          /* Tvaraff[ij]=  4    3                                            1  */
          /* Tmodelind[ij]=2    3                                            9  */
          /* TmodelInvind[ij]=2 1                                            1  */
        Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/         Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/
        Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */         Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */
        TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */         TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */
Line 6063  void  concatwav(int wav[], int **dh, int Line 6140  void  concatwav(int wav[], int **dh, int
    } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */     } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */
    /* ij--; */     /* ij--; */
    /* cptcoveff=ij; /\*Number of total covariates*\/ */     /* cptcoveff=ij; /\*Number of total covariates*\/ */
    *cptcov=ij; /* cptcov= Number of total real effective covariates: effective (used as cptcoveff in other functions)     *cptcov=ij; /* cptcov= Number of total real effective simple dummies (fixed or time  arying) effective (used as cptcoveff in other functions)
                 * because they can be excluded from the model and real                  * because they can be excluded from the model and real
                 * if in the model but excluded because missing values, but how to get k from ij?*/                  * if in the model but excluded because missing values, but how to get k from ij?*/
    for(j=ij+1; j<= cptcovt; j++){     for(j=ij+1; j<= cptcovt; j++){
Line 6196  void  concatwav(int wav[], int **dh, int Line 6273  void  concatwav(int wav[], int **dh, int
   /* Covariances of health expectancies eij and of total life expectancies according    /* Covariances of health expectancies eij and of total life expectancies according
      to initial status i, ei. .       to initial status i, ei. .
   */    */
     /* Very time consuming function, but already optimized with precov */
   int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji;    int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji;
   int nhstepma, nstepma; /* Decreasing with age */    int nhstepma, nstepma; /* Decreasing with age */
   double age, agelim, hf;    double age, agelim, hf;
Line 6464  void  concatwav(int wav[], int **dh, int Line 6542  void  concatwav(int wav[], int **dh, int
    pstamp(ficresprobmorprev);     pstamp(ficresprobmorprev);
    fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);     fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
    fprintf(ficresprobmorprev,"# Selected quantitative variables and dummies");     fprintf(ficresprobmorprev,"# Selected quantitative variables and dummies");
    for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */     for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ /* To be done*/
      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<=cptcoveff;j++)      for(j=1;j<=cptcoveff;j++) 
      fprintf(ficresprobmorprev,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(ij,TnsdVar[Tvaraff[j]])]);       fprintf(ficresprobmorprev," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(ij,TnsdVar[Tvaraff[j]])]);
    fprintf(ficresprobmorprev,"\n");     fprintf(ficresprobmorprev,"\n");
   
    fprintf(ficresprobmorprev,"# Age cov=%-d",ij);     fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
Line 7099  To be simple, these graphs help to under Line 7177  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 \n",  TKresult[nres], j1, nres, cptcovn, cptcoveff, tj);       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 7107  To be simple, these graphs help to under Line 7185  To be simple, these graphs help to under
      /* for(nres=1;nres <=1; nres++){ /\* For each resultline *\/ */       /* for(nres=1;nres <=1; nres++){ /\* For each resultline *\/ */
      /* /\* for(nres=1;nres <=nresult; nres++){ /\\* For each resultline *\\/ *\/ */       /* /\* for(nres=1;nres <=nresult; nres++){ /\\* For each resultline *\\/ *\/ */
      if  (cptcovn>0) {       if  (cptcovn>0) {
        fprintf(ficresprob, "\n#********** Variable ");          fprintf(ficresprob, "\n#********** Variable ");
        for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);  
        fprintf(ficresprob, "**********\n#\n");  
        fprintf(ficresprobcov, "\n#********** Variable ");          fprintf(ficresprobcov, "\n#********** Variable "); 
        for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);         fprintf(ficgp, "\n#********** Variable ");
          fprintf(fichtmcov, "\n<hr  size=\"2\" color=\"#EC5E5E\">********** Variable "); 
          fprintf(ficresprobcor, "\n#********** Variable ");    
   
          /* Including quantitative variables of the resultline to be done */
          for (z1=1; z1<=cptcovs; z1++){ /* Loop on each variable of this resultline  */
            printf("Varprob modelresult[%d][%d]=%d model=%s \n",nres, z1, modelresult[nres][z1], model);
            fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=%s \n",nres, z1, modelresult[nres][z1], model);
            /* fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=%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(Fixed[modelresult[nres][z1]]==0){ /* Fixed referenced to model equation */
                fprintf(ficresprob,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
                fprintf(ficresprobcov,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
                fprintf(ficgp,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
                fprintf(fichtmcov,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
                fprintf(ficresprobcor,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
                fprintf(ficresprob,"fixed ");
                fprintf(ficresprobcov,"fixed ");
                fprintf(ficgp,"fixed ");
                fprintf(fichtmcov,"fixed ");
                fprintf(ficresprobcor,"fixed ");
              }else{
                fprintf(ficresprob,"varyi ");
                fprintf(ficresprobcov,"varyi ");
                fprintf(ficgp,"varyi ");
                fprintf(fichtmcov,"varyi ");
                fprintf(ficresprobcor,"varyi ");
              }
            }else if(Dummy[modelresult[nres][z1]]==1){ /* Quanti variable */
              /* For each selected (single) quantitative value */
              fprintf(ficresprob," V%d=%f ",Tvqresult[nres][z1],Tqresult[nres][z1]);
              if(Fixed[modelresult[nres][z1]]==0){ /* Fixed */
                fprintf(ficresprob,"fixed ");
                fprintf(ficresprobcov,"fixed ");
                fprintf(ficgp,"fixed ");
                fprintf(fichtmcov,"fixed ");
                fprintf(ficresprobcor,"fixed ");
              }else{
                fprintf(ficresprob,"varyi ");
                fprintf(ficresprobcov,"varyi ");
                fprintf(ficgp,"varyi ");
                fprintf(fichtmcov,"varyi ");
                fprintf(ficresprobcor,"varyi ");
              }
            }else{
              printf("Error in varprob() Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=V%d cptcovs=%d, cptcoveff=%d \n", nres, z1, Dummy[modelresult[nres][z1]],nres,z1,modelresult[nres][z1],cptcovs, cptcoveff);  /* end if dummy  or quanti */
              fprintf(ficlog,"Error in varprob() Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=V%d cptcovs=%d, cptcoveff=%d \n", nres, z1, Dummy[modelresult[nres][z1]],nres,z1,modelresult[nres][z1],cptcovs, cptcoveff);  /* end if dummy  or quanti */
              exit(1);
            }
          } /* End loop on variable of this resultline */
          /* for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); */
          fprintf(ficresprob, "**********\n#\n");
        fprintf(ficresprobcov, "**********\n#\n");         fprintf(ficresprobcov, "**********\n#\n");
                           
        fprintf(ficgp, "\n#********** Variable ");   
        for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);  
        fprintf(ficgp, "**********\n#\n");         fprintf(ficgp, "**********\n#\n");
                           
                           
        fprintf(fichtmcov, "\n<hr  size=\"2\" color=\"#EC5E5E\">********** Variable ");   
        /* for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); */  
        for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtmcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);  
        fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");
                           
        fprintf(ficresprobcor, "\n#********** Variable ");      
        for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);  
        fprintf(ficresprobcor, "**********\n#");             fprintf(ficresprobcor, "**********\n#");    
        if(invalidvarcomb[j1]){         if(invalidvarcomb[j1]){
          fprintf(ficgp,"\n#Combination (%d) ignored because no cases \n",j1);            fprintf(ficgp,"\n#Combination (%d) ignored because no cases \n",j1); 
Line 7137  To be simple, these graphs help to under Line 7253  To be simple, these graphs help to under
      trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);       trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
      gp=vector(1,(nlstate)*(nlstate+ndeath));       gp=vector(1,(nlstate)*(nlstate+ndeath));
      gm=vector(1,(nlstate)*(nlstate+ndeath));       gm=vector(1,(nlstate)*(nlstate+ndeath));
      for (age=bage; age<=fage; age ++){        for (age=bage; age<=fage; age ++){ /* Fo each age we feed the model equation with covariates, using precov as in hpxij() ? */
        cov[2]=age;         cov[2]=age;
        if(nagesqr==1)         if(nagesqr==1)
          cov[3]= age*age;           cov[3]= age*age;
        /* for (k=1; k<=cptcovn;k++) { */         /* New code end of combination but for each resultline */
        /*        cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; */         for(k1=1;k1<=cptcovt;k1++){ /* loop on model equation (including products) */ 
        for (k=1; k<=nsd;k++) { /* For single dummy covariates only */           if(Typevar[k1]==1){ /* A product with age */
          /* Here comes the value of the covariate 'j1' after renumbering k with single dummy covariates */             cov[2+nagesqr+k1]=precov[nres][k1]*cov[2];
          cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(j1,TnsdVar[TvarsD[k]])];  
          /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4  
                                                                     * 1  1 1 1 1  
                                                                     * 2  2 1 1 1  
                                                                     * 3  1 2 1 1  
                                                                     */  
          /* nbcode[1][1]=0 nbcode[1][2]=1;*/  
        }  
        /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */  
        /* ) p nbcode[Tvar[Tage[k]]][(1 & (ij-1) >> (k-1))+1] */  
        /*for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */  
        for (k=1; k<=cptcovage;k++){  /* For product with age */  
          if(Dummy[Tage[k]]==2){ /* dummy with age */  
            cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(j1,TnsdVar[Tvar[Tage[k]]])]*cov[2];  
            /* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */  
          } else if(Dummy[Tage[k]]==3){ /* quantitative with age */  
            printf("Internal IMaCh error, don't know which value for quantitative covariate with age, Tage[k]%d, k=%d, Tvar[Tage[k]]=V%d, age=%d\n",Tage[k],k ,Tvar[Tage[k]], (int)cov[2]);  
            /* cov[2+nagesqr+Tage[k]]=meanq[k]/idq[k]*cov[2];/\* Using the mean of quantitative variable Tvar[Tage[k]] /\* Tqresult[nres][k]; *\/ */  
            /* exit(1); */  
            /* cov[++k1]=Tqresult[nres][k];  */  
          }  
          /* cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */  
        }  
        for (k=1; k<=cptcovprod;k++){/* For product without age */  
          if(Dummy[Tvard[k][1]]==0){  
            if(Dummy[Tvard[k][2]]==0){  
              cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,TnsdVar[Tvard[k][1]])] * nbcode[Tvard[k][2]][codtabm(j1,TnsdVar[Tvard[k][2]])];  
              /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */  
            }else{ /* Should we use the mean of the quantitative variables? */  
              cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,TnsdVar[Tvard[k][1]])] * Tqresult[nres][resultmodel[nres][k]];  
              /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; */  
            }  
          }else{           }else{
            if(Dummy[Tvard[k][2]]==0){             cov[2+nagesqr+k1]=precov[nres][k1];
              cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(j1,TnsdVar[Tvard[k][2]])] * Tqinvresult[nres][TnsdVar[Tvard[k][1]]];  
              /* cov[++k1]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; */  
            }else{  
              cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][TnsdVar[Tvard[k][1]]]*  Tqinvresult[nres][TnsdVar[Tvard[k][2]]];  
              /* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]]; */  
            }  
          }           }
          /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; */         }/* End of loop on model equation */
        }                          /* Old code */
          /* /\* for (k=1; k<=cptcovn;k++) { *\/ */
          /* /\*    cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; *\/ */
          /* for (k=1; k<=nsd;k++) { /\* For single dummy covariates only *\/ */
          /*        /\* Here comes the value of the covariate 'j1' after renumbering k with single dummy covariates *\/ */
          /*        cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(j1,TnsdVar[TvarsD[k]])]; */
          /*        /\*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*\//\* j1 1 2 3 4 */
          /*                                                                   * 1  1 1 1 1 */
          /*                                                                   * 2  2 1 1 1 */
          /*                                                                   * 3  1 2 1 1 */
          /*                                                                   *\/ */
          /*        /\* nbcode[1][1]=0 nbcode[1][2]=1;*\/ */
          /* } */
          /* /\* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 *\/ */
          /* /\* ) p nbcode[Tvar[Tage[k]]][(1 & (ij-1) >> (k-1))+1] *\/ */
          /* /\*for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; *\/ */
          /* for (k=1; k<=cptcovage;k++){  /\* For product with age *\/ */
          /*        if(Dummy[Tage[k]]==2){ /\* dummy with age *\/ */
          /*          cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(j1,TnsdVar[Tvar[Tage[k]]])]*cov[2]; */
          /*          /\* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; *\/ */
          /*        } else if(Dummy[Tage[k]]==3){ /\* quantitative with age *\/ */
          /*          printf("Internal IMaCh error, don't know which value for quantitative covariate with age, Tage[k]%d, k=%d, Tvar[Tage[k]]=V%d, age=%d\n",Tage[k],k ,Tvar[Tage[k]], (int)cov[2]); */
          /*          /\* cov[2+nagesqr+Tage[k]]=meanq[k]/idq[k]*cov[2];/\\* Using the mean of quantitative variable Tvar[Tage[k]] /\\* Tqresult[nres][k]; *\\/ *\/ */
          /*          /\* exit(1); *\/ */
          /*          /\* cov[++k1]=Tqresult[nres][k];  *\/ */
          /*        } */
          /*        /\* cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; *\/ */
          /* } */
          /* for (k=1; k<=cptcovprod;k++){/\* For product without age *\/ */
          /*        if(Dummy[Tvard[k][1]]==0){ */
          /*          if(Dummy[Tvard[k][2]]==0){ */
          /*            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,TnsdVar[Tvard[k][1]])] * nbcode[Tvard[k][2]][codtabm(j1,TnsdVar[Tvard[k][2]])]; */
          /*            /\* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; *\/ */
          /*          }else{ /\* Should we use the mean of the quantitative variables? *\/ */
          /*            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,TnsdVar[Tvard[k][1]])] * Tqresult[nres][resultmodel[nres][k]]; */
          /*            /\* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; *\/ */
          /*          } */
          /*        }else{ */
          /*          if(Dummy[Tvard[k][2]]==0){ */
          /*            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(j1,TnsdVar[Tvard[k][2]])] * Tqinvresult[nres][TnsdVar[Tvard[k][1]]]; */
          /*            /\* cov[++k1]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; *\/ */
          /*          }else{ */
          /*            cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][TnsdVar[Tvard[k][1]]]*  Tqinvresult[nres][TnsdVar[Tvard[k][2]]]; */
          /*            /\* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]]; *\/ */
          /*          } */
          /*        } */
          /*        /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; *\/ */
          /* } */                  
 /* For each age and combination of dummy covariates we slightly move the parameters of delti in order to get the gradient*/                       /* For each age and combination of dummy covariates we slightly move the parameters of delti in order to get the gradient*/                     
        for(theta=1; theta <=npar; theta++){         for(theta=1; theta <=npar; theta++){
          for(i=1; i<=npar; i++)           for(i=1; i<=npar; i++)
Line 10265  int decoderesult( char resultline[], int Line 10390  int decoderesult( char resultline[], int
   int j=0, k=0, k1=0, k2=0, k3=0, k4=0, match=0, k2q=0, k3q=0, k4q=0;    int j=0, k=0, k1=0, k2=0, k3=0, k4=0, match=0, k2q=0, k3q=0, k4q=0;
   char resultsav[MAXLINE];    char resultsav[MAXLINE];
   /* int resultmodel[MAXLINE]; */    /* int resultmodel[MAXLINE]; */
   int modelresult[MAXLINE];    /* int modelresult[MAXLINE]; */
   char stra[80], strb[80], strc[80], strd[80],stre[80];    char stra[80], strb[80], strc[80], strd[80],stre[80];
   
   removefirstspace(&resultline);    removefirstspace(&resultline);
Line 10274  int decoderesult( char resultline[], int Line 10399  int decoderesult( char resultline[], int
   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'=' */      j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' in this resultline */
   }    }
   if(j == 0){ /* Resultline but no = */    if(j == 0){ /* Resultline but no = */
     TKresult[nres]=0; /* Combination for the nresult and the model */      TKresult[nres]=0; /* Combination for the nresult and the model */
     return (0);      return (0);
   }    }
   if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */    if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */
     printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of variables used in the model line, %s.\n",j, cptcovs, model);      printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model);
     fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of variables used in the model line, %s.\n",j, cptcovs, model);      fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model);
     /* return 1;*/      /* return 1;*/
   }    }
   for(k=1; k<=j;k++){ /* Loop on any covariate of the result line */    for(k=1; k<=j;k++){ /* Loop on any covariate of the RESULT LINE */
     if(nbocc(resultsav,'=') >1){      if(nbocc(resultsav,'=') >1){
       cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' ' (stra is the rest of the resultline to be analyzed in the next loop *//*     resultsav= "V4=1 V5=25.1 V3=0" stra= "V5=25.1 V3=0" strb= "V4=1" */        cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' ' (stra is the rest of the resultline to be analyzed in the next loop *//*     resultsav= "V4=1 V5=25.1 V3=0" stra= "V5=25.1 V3=0" strb= "V4=1" */
       /* If resultsav= "V4= 1 V5=25.1 V3=0" with a blank then strb="V4=" and stra="1 V5=25.1 V3=0" */        /* If resultsav= "V4= 1 V5=25.1 V3=0" with a blank then strb="V4=" and stra="1 V5=25.1 V3=0" */
Line 10309  int decoderesult( char resultline[], int Line 10434  int decoderesult( char resultline[], int
   }    }
   /* Checking for missing or useless values in comparison of current model needs */    /* Checking for missing or useless values in comparison of current model needs */
   /* Feeds resultmodel[nres][k1]=k2 for k1th product covariate with age in the model equation fed by the index k2 of the resutline*/    /* Feeds resultmodel[nres][k1]=k2 for k1th product covariate with age in the model equation fed by the index k2 of the resutline*/
   for(k1=1; k1<= cptcovt ;k1++){ /* Loop on model. model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */    for(k1=1; k1<= cptcovt ;k1++){ /* Loop on MODEL LINE V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
     if(Typevar[k1]==0){ /* Single covariate in model */      if(Typevar[k1]==0){ /* Single covariate in model */
       /* 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product */        /* 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product */
       match=0;        match=0;
       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(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[k2]=k1;/* modelresult[2]=1 modelresult[1]=2  modelresult[3]=3  modelresult[6]=4 modelresult[9]=5 */            modelresult[nres][k2]=k1;/* modelresult[2]=1 modelresult[1]=2  modelresult[3]=3  modelresult[6]=4 modelresult[9]=5 */
           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 10330  int decoderesult( char resultline[], int Line 10455  int decoderesult( char resultline[], int
       match=0;        match=0;
       for(k2=1; k2 <=j;k2++){/* Loop on resultline.  jth occurence of = signs in the result line. In result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */        for(k2=1; k2 <=j;k2++){/* Loop on resultline.  jth occurence of = signs in the result line. In result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
         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[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 */
Line 10376  int decoderesult( char resultline[], int Line 10501  int decoderesult( char resultline[], int
   }/* End loop cptcovt */    }/* End loop cptcovt */
   /* Checking for missing or useless values in comparison of current model needs */    /* Checking for missing or useless values in comparison of current model needs */
   /* Feeds resultmodel[nres][k1]=k2 for single covariate (k1) in the model equation */    /* Feeds resultmodel[nres][k1]=k2 for single covariate (k1) in the model equation */
   for(k2=1; k2 <=j;k2++){ /* Loop on resultline variables: result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */    for(k2=1; k2 <=j;k2++){ /* j or cptcovs is the number of single covariates used either in the model line as well as in the result line (dummy or quantitative)
                              * Loop on resultline variables: result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
     match=0;      match=0;
     for(k1=1; k1<= cptcovt ;k1++){ /* loop on model: model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */      for(k1=1; k1<= cptcovt ;k1++){ /* loop on model: model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
       if(Typevar[k1]==0){ /* Single only */        if(Typevar[k1]==0){ /* Single only */
         if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=4 == Tvarsel[1]=4   */          if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=4 == Tvarsel[1]=4   */
           resultmodel[nres][k1]=k2;  /* k1th position in the model equation corresponds to k2th position in the result line. resultmodel[2]=1 resultmodel[1]=2  resultmodel[3]=3  resultmodel[6]=4 resultmodel[9]=5 */            resultmodel[nres][k1]=k2;  /* k1th position in the model equation corresponds to k2th position in the result line. resultmodel[2]=1 resultmodel[1]=2  resultmodel[3]=3  resultmodel[6]=4 resultmodel[9]=5 */
             modelresult[nres][k2]=k1; /* k1th position in the model equation corresponds to k2th position in the result line. modelresult[1]=2 modelresult[2]=1  modelresult[3]=3  remodelresult[4]=6 modelresult[5]=9 */
           ++match;            ++match;
         }          }
       }        }
Line 10396  int decoderesult( char resultline[], int Line 10523  int decoderesult( char resultline[], int
       return 1;        return 1;
     }      }
   }    }
           /* cptcovres=j /\* Number of variables in the resultline is equal to cptcovs and thus useless *\/     */
   /* We need to deduce which combination number is chosen and save quantitative values */    /* We need to deduce which combination number is chosen and save quantitative values */
   /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */    /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
   /* nres=1st result line: V4=1 V5=25.1 V3=0  V2=8 V1=1 */    /* nres=1st result line: V4=1 V5=25.1 V3=0  V2=8 V1=1 */
Line 10415  int decoderesult( char resultline[], int Line 10542  int decoderesult( char resultline[], int
   /* V(Tvqresult)=Tqresult V5=25.1 V2=8 Tqresult[nres=1][1]=25.1 */    /* V(Tvqresult)=Tqresult V5=25.1 V2=8 Tqresult[nres=1][1]=25.1 */
   /* V5*age V5 known which value for nres?  */    /* V5*age V5 known which value for nres?  */
   /* Tqinvresult[2]=8 Tqinvresult[1]=25.1  */    /* Tqinvresult[2]=8 Tqinvresult[1]=25.1  */
   for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* loop k1 on position in the model line (excluding product) */    for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* cptcovt number of covariates (excluding 1 and age or age*age) in the MODEL equation.
                                                      * loop on position k1 in the MODEL LINE */
     /* 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 */      if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Dummy and Single, k1 is sorting according to MODEL, but k3 to resultline */
        /* k4+1= 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 */
       /* resultmodel[nres][k1]=k3: k1th position in the model correspond to the k3 position in the resultline                        */        /* resultmodel[nres][k1]=k3: k1th position in the model correspond to the k3 position in the resultline                        */
Line 10428  int decoderesult( char resultline[], int Line 10556  int decoderesult( char resultline[], int
       /* Tvarsel[k3]: Name of the variable at the k3th position in the result line.                                                  */        /* Tvarsel[k3]: Name of the variable at the k3th position in the result line.                                                  */
       /* Tvalsel[k3]: Value of the variable at the k3th position in the result line.                                                 */        /* Tvalsel[k3]: Value of the variable at the k3th position in the result line.                                                 */
       /* Tresult[nres][result_position]= value of the dummy variable at the result_position in the nres resultline                   */        /* Tresult[nres][result_position]= value of the dummy variable at the result_position in the nres resultline                   */
       /* Tvresult[nres][result_position]= id of the dummy variable at the result_position in the nres resultline                     */        /* Tvresult[nres][result_position]= name of the dummy variable at the result_position in the nres resultline                     */
       /* Tinvresult[nres][Name of a dummy variable]= value of the variable in the result line                                        */        /* Tinvresult[nres][Name of a dummy variable]= value of the variable in the result line                                        */
       /* TinvDoQresult[nres][Name of a Dummy or Q variable]= value of the variable in the result line                                                      */        /* TinvDoQresult[nres][Name of a Dummy or Q variable]= value of the variable in the result line                                                      */
       k3= resultmodel[nres][k1]; /* From position k1 in model get position k3 in result line */        k3= resultmodel[nres][k1]; /* From position k1 in model get position k3 in result line */
       /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/        /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/
       k2=(int)Tvarsel[k3]; /* from position k3 in resultline get name k2: 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]; /* from position k3 in resultline get name k2: nres=1 k1=2=>k3=1 Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 (V4); k1=3=>k3=2 Tvarsel[2]=3 (V3)*/
       k+=Tvalsel[k3]*pow(2,k4);  /* nres=1 k1=2 Tvalsel[1]=1 (V4=1); k1=3 k3=2 Tvalsel[2]=0 (V3=0) */        k+=Tvalsel[k3]*pow(2,k4);  /* nres=1 k1=2 Tvalsel[1]=1 (V4=1); k1=3 k3=2 Tvalsel[2]=0 (V3=0) */
       TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Stores the value into the name of the variable. */        TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* TinvDoQresult[nres][Name]=Value; stores the value into the name of the variable. */
       /* Tinvresult[nres][4]=1 */        /* Tinvresult[nres][4]=1 */
       Tresult[nres][k4+1]=Tvalsel[k3];/* Tresult[nres=2][1]=1(V4=1)  Tresult[nres=2][2]=0(V3=0) */        /* Tresult[nres][k4+1]=Tvalsel[k3];/\* Tresult[nres=2][1]=1(V4=1)  Tresult[nres=2][2]=0(V3=0) *\/ */
       Tvresult[nres][k4+1]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */        Tresult[nres][k3]=Tvalsel[k3];/* Tresult[nres=2][1]=1(V4=1)  Tresult[nres=2][2]=0(V3=0) */
         /* Tvresult[nres][k4+1]=(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];        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 */
Line 10450  int decoderesult( char resultline[], int Line 10580  int decoderesult( char resultline[], int
       k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 5 =k3q */        k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 5 =k3q */
       k2q=(int)Tvarsel[k3q]; /*  Name of variable at k3q th position in the resultline */        k2q=(int)Tvarsel[k3q]; /*  Name of variable at k3q th position in the resultline */
       /* Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */        /* Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */
       Tqresult[nres][k4q+1]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */        /* Tqresult[nres][k4q+1]=Tvalsel[k3q]; /\* Tqresult[nres][1]=25.1 *\/ */
       Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */        /* Tvresult[nres][k4q+1]=(int)Tvarsel[k3q];/\* Tvresult[nres][1]=4 Tvresult[nres][3]=1 *\/ */
         /* Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /\* Tvqresult[nres][1]=5 *\/ */
         Tqresult[nres][k3q]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */
         Tvresult[nres][k3q]=(int)Tvarsel[k3q];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */
         Tvqresult[nres][k3q]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */
       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];
Line 10463  int decoderesult( char resultline[], int Line 10597  int decoderesult( char resultline[], int
               
       k3= resultmodel[nres][k1]; /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/        k3= resultmodel[nres][k1]; /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/
       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]; /* Tinvresult[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]; /* Tqinvresult[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]]);
Line 10481  int decoderesult( char resultline[], int Line 10615  int decoderesult( char resultline[], int
     }      }
   }    }
       
   TKresult[nres]=++k; /* Combination for the nresult and the model */    TKresult[nres]=++k; /* Number of combinations of dummies for the nresult and the model =Tvalsel[k3]*pow(2,k4) + 1*/
   return (0);    return (0);
 }  }
   
Line 10641  int decodemodel( char model[], int lasto Line 10775  int decodemodel( char model[], int lasto
             Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but              Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but
                                                 because this model-covariate is a construction we invent a new column                                                  because this model-covariate is a construction we invent a new column
                                                 which is after existing variables ncovcol+nqv+ntv+nqtv + k1                                                  which is after existing variables ncovcol+nqv+ntv+nqtv + k1
                                                 If already ncovcol=4 and model=V2 + V1 +V1*V4 +age*V3 +V3*V2                                                  If already ncovcol=4 and model= V2 + V1 + V1*V4 + age*V3 + V3*V2
                                                 thus after V4 we invent V5 and V6 because age*V3 will be computed in 4                                                  thus after V4 we invent V5 and V6 because age*V3 will be computed in 4
                                                 Tvar[3=V1*V4]=4+1=5 Tvar[5=V3*V2]=4 + 2= 6, Tvar[4=age*V3]=4 etc */                                                  Tvar[3=V1*V4]=4+1=5 Tvar[5=V3*V2]=4 + 2= 6, Tvar[4=age*V3]=4 etc */
               /* Please remark that the new variables are model dependent */
               /* If we have 4 variable but the model uses only 3, like in
                * model= V1 + age*V1 + V2 + V3 + age*V2 + age*V3 + V1*V2 + V1*V3
                *  k=     1     2       3   4     5        6        7       8
                * Tvar[k]=1     1       2   3     2        3       (5       6) (and not 4 5 because of V4 missing)
                * Tage[kk]    [1]= 2           [2]=5      [3]=6                  kk=1 to cptcovage=3
                * Tvar[Tage[kk]][1]=2          [2]=2      [3]=3
                */
             Typevar[k]=2;  /* 2 for double fixed dummy covariates */              Typevar[k]=2;  /* 2 for double fixed dummy covariates */
             cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */              cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
             Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */              Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */
Line 10746  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 10888  Dummy[k] 0=dummy (0 1), 1 quantitative (
       modell[k].maintype= FTYPE;        modell[k].maintype= FTYPE;
       modell[k].subtype= FQ;        modell[k].subtype= FQ;
       nsq++;        nsq++;
       TvarsQ[nsq]=Tvar[k];        TvarsQ[nsq]=Tvar[k]; /* Gives the variable name (extended to products) of first nsq simple quantitative covariates (fixed or time vary see below */
       TvarsQind[nsq]=k;        TvarsQind[nsq]=k;    /* Gives the position in the model equation of the first nsq simple quantitative covariates (fixed or time vary) */
       ncovf++;        ncovf++;
       TvarF[ncovf]=Tvar[k];        TvarF[ncovf]=Tvar[k];
       TvarFind[ncovf]=k;        TvarFind[ncovf]=k;
Line 10778  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 10920  Dummy[k] 0=dummy (0 1), 1 quantitative (
       modell[k].subtype= VQ;        modell[k].subtype= VQ;
       ncovv++; /* Only simple time varying variables */        ncovv++; /* Only simple time varying variables */
       nsq++;        nsq++;
       TvarsQ[nsq]=Tvar[k]; /* k=1 Tvar=5 nsq=1 TvarsQ[1]=5 */        TvarsQ[nsq]=Tvar[k]; /* k=1 Tvar=5 nsq=1 TvarsQ[1]=5 */ /* Gives the variable name (extended to products) of first nsq simple quantitative covariates (fixed or time vary here) */
       TvarsQind[nsq]=k; /* For single quantitative covariate gives the model position of each single quantitative covariate */        TvarsQind[nsq]=k; /* For single quantitative covariate gives the model position of each single quantitative covariate *//* Gives the position in the model equation of the first nsq simple quantitative covariates (fixed or time vary) */
       TvarV[ncovv]=Tvar[k];        TvarV[ncovv]=Tvar[k];
       TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */        TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */
       TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */        TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */
Line 11506  int back_prevalence_limit(double *p, dou Line 11648  int back_prevalence_limit(double *p, dou
     
 int hPijx(double *p, int bage, int fage){  int hPijx(double *p, int bage, int fage){
     /*------------- h Pij x at various ages ------------*/      /*------------- h Pij x at various ages ------------*/
     /* to be optimized with precov */
   int stepsize;    int stepsize;
   int agelim;    int agelim;
   int hstepm;    int hstepm;
Line 11583  int hPijx(double *p, int bage, int fage) Line 11725  int hPijx(double *p, int bage, int fage)
     
  int hBijx(double *p, int bage, int fage, double ***prevacurrent){   int hBijx(double *p, int bage, int fage, double ***prevacurrent){
     /*------------- h Bij x at various ages ------------*/      /*------------- h Bij x at various ages ------------*/
       /* To be optimized with precov */
   int stepsize;    int stepsize;
   /* int agelim; */    /* int agelim; */
         int ageminl;          int ageminl;
Line 11722  int main(int argc, char *argv[]) Line 11864  int main(int argc, char *argv[])
       
   char pathr[MAXLINE], pathimach[MAXLINE];     char pathr[MAXLINE], pathimach[MAXLINE]; 
   char *tok, *val; /* pathtot */    char *tok, *val; /* pathtot */
   int firstobs=1, lastobs=10; /* nobs = lastobs-firstobs declared globally ;*/    /* int firstobs=1, lastobs=10; /\* nobs = lastobs-firstobs declared globally ;*\/ */
   int c,  h , cpt, c2;    int c,  h , cpt, c2;
   int jl=0;    int jl=0;
   int i1, j1, jk, stepsize=0;    int i1, j1, jk, stepsize=0;
Line 12345  Please run with mle=-1 to get a correct Line 12487  Please run with mle=-1 to get a correct
   mint=matrix(1,maxwav,firstobs,lastobs);    mint=matrix(1,maxwav,firstobs,lastobs);
   anint=matrix(1,maxwav,firstobs,lastobs);    anint=matrix(1,maxwav,firstobs,lastobs);
   s=imatrix(1,maxwav+1,firstobs,lastobs); /* s[i][j] health state for wave i and individual j */    s=imatrix(1,maxwav+1,firstobs,lastobs); /* s[i][j] health state for wave i and individual j */
   printf("BUG ncovmodel=%d NCOVMAX=%d 2**ncovmodel=%f BUG\n",ncovmodel,NCOVMAX,pow(2,ncovmodel));    /* printf("BUG ncovmodel=%d NCOVMAX=%d 2**ncovmodel=%f BUG\n",ncovmodel,NCOVMAX,pow(2,ncovmodel)); */
   tab=ivector(1,NCOVMAX);    tab=ivector(1,NCOVMAX);
   ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */    ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */
   ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */    ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */
Line 12365  Please run with mle=-1 to get a correct Line 12507  Please run with mle=-1 to get a correct
   Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */    Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */
   TvarsDind=ivector(1,NCOVMAX); /*  */    TvarsDind=ivector(1,NCOVMAX); /*  */
   TnsdVar=ivector(1,NCOVMAX); /*  */    TnsdVar=ivector(1,NCOVMAX); /*  */
       /* for(i=1; i<=NCOVMAX;i++) TnsdVar[i]=3; */
   TvarsD=ivector(1,NCOVMAX); /*  */    TvarsD=ivector(1,NCOVMAX); /*  */
   TvarsQind=ivector(1,NCOVMAX); /*  */    TvarsQind=ivector(1,NCOVMAX); /*  */
   TvarsQ=ivector(1,NCOVMAX); /*  */    TvarsQ=ivector(1,NCOVMAX); /*  */
Line 12487  Please run with mle=-1 to get a correct Line 12630  Please run with mle=-1 to get a correct
   Ndum =ivector(-1,NCOVMAX);      Ndum =ivector(-1,NCOVMAX);  
   cptcoveff=0;    cptcoveff=0;
   if (ncovmodel-nagesqr > 2 ){ /* That is if covariate other than cst, age and age*age */    if (ncovmodel-nagesqr > 2 ){ /* That is if covariate other than cst, age and age*age */
     tricode(&cptcoveff,Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */      tricode(&cptcoveff,Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; as well as calculate cptcoveff or number of total effective dummy covariates*/
   }    }
       
   ncovcombmax=pow(2,cptcoveff);    ncovcombmax=pow(2,cptcoveff);
Line 12626  Title=%s <br>Datafile=%s Firstpass=%d La Line 12769  Title=%s <br>Datafile=%s Firstpass=%d La
           optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);            optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }    }
   
   fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C)  2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-2013-2022-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br>  \    fprintf(fichtm,"<html><head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n\
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <title>IMaCh %s</title></head>\n\
    <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n\
   <font size=\"3\">Sponsored by Copyright (C)  2002-2015 <a href=http://www.ined.fr>INED</a>\
   -EUROREVES-Institut de longévité-2013-2022-Japan Society for the Promotion of Sciences 日本学術振興会 \
   (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - \
   <a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br> \n", optionfilehtm);
     
     fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\"> \n\
 <font size=\"2\">IMaCh-%s <br> %s</font> \  <font size=\"2\">IMaCh-%s <br> %s</font> \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n\  This file: <a href=\"%s\">%s</a>Title=%s <br>Datafile=<a href=\"%s\">%s</a> Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n\
 \n\  \n\
 <hr  size=\"2\" color=\"#EC5E5E\">\  <hr  size=\"2\" color=\"#EC5E5E\">\
  <ul><li><h4>Parameter files</h4>\n\   <ul><li><h4>Parameter files</h4>\n\
Line 12639  Title=%s <br>Datafile=%s Firstpass=%d La Line 12789  Title=%s <br>Datafile=%s Firstpass=%d La
  - Log file of the run: <a href=\"%s\">%s</a><br>\n\   - Log file of the run: <a href=\"%s\">%s</a><br>\n\
  - Gnuplot file name: <a href=\"%s\">%s</a><br>\n\   - Gnuplot file name: <a href=\"%s\">%s</a><br>\n\
  - Date and time at start: %s</ul>\n",\   - Date and time at start: %s</ul>\n",\
           optionfilehtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\            version,fullversion,optionfilehtm,optionfilehtm,title,datafile,datafile,firstpass,lastpass,stepm, weightopt, model, \
           optionfilefiname,optionfilext,optionfilefiname,optionfilext,\            optionfilefiname,optionfilext,optionfilefiname,optionfilext,\
           fileres,fileres,\            fileres,fileres,\
           filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);            filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);
Line 12953  Please run with mle=-1 to get a correct Line 13103  Please run with mle=-1 to get a correct
     globpr=1; /* again, to print the individual contributions using computed gpimx and gsw */      globpr=1; /* again, to print the individual contributions using computed gpimx and gsw */
     likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */      likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
     printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);      printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
             /* exit(0); */
     for (k=1; k<=npar;k++)      for (k=1; k<=npar;k++)
       printf(" %d %8.5f",k,p[k]);        printf(" %d %8.5f",k,p[k]);
     printf("\n");      printf("\n");
Line 13631  Please run with mle=-1 to get a correct Line 13782  Please run with mle=-1 to get a correct
         printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);          printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
       }        }
       for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */        for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
         printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);          printf(" V%d=%f ",TvarsQ[j], TinvDoQresult[nres][TvarsQ[j]]); /* TvarsQ[j] gives the name of the jth quantitative (fixed or time v) */
         fprintf(ficreseij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);          fprintf(ficreseij,"V%d=%f ",TvarsQ[j], TinvDoQresult[nres][TvarsQ[j]]);
       }        }
       fprintf(ficreseij,"******\n");        fprintf(ficreseij,"******\n");
       printf("******\n");        printf("******\n");
Line 13650  Please run with mle=-1 to get a correct Line 13801  Please run with mle=-1 to get a correct
   
                                   
     /*---------- State-specific expectancies and variances ------------*/      /*---------- State-specific expectancies and variances ------------*/
                       /* Should be moved in a function */         
     strcpy(filerest,"T_");      strcpy(filerest,"T_");
     strcat(filerest,fileresu);      strcat(filerest,fileresu);
     if((ficrest=fopen(filerest,"w"))==NULL) {      if((ficrest=fopen(filerest,"w"))==NULL) {
Line 13689  Please run with mle=-1 to get a correct Line 13840  Please run with mle=-1 to get a correct
     i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */      i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
     if (cptcovn < 1){i1=1;}      if (cptcovn < 1){i1=1;}
           
     for(nres=1; nres <= nresult; nres++) /* For each resultline */      for(nres=1; nres <= nresult; nres++) /* For each resultline, find the combination and output results according to the values of dummies and then quanti.  */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */      for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying. For each nres and each value at position k
       if(i1 != 1 && TKresult[nres]!= k)                            * we know Tresult[nres][result_position]= value of the dummy variable at the result_position in the nres resultline
                             * Tvqresult[nres][result_position]= id of the variable at the result_position in the nres resultline 
                             * and Tqresult[nres][result_position]= value of the variable at the result_position in the nres resultline */
         /* */
         if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */
         continue;          continue;
       printf("\n# model %s \n#****** Result for:", model);        printf("\n# model %s \n#****** Result for:", model);
       fprintf(ficrest,"\n# model %s \n#****** Result for:", model);        fprintf(ficrest,"\n# model %s \n#****** Result for:", model);
       fprintf(ficlog,"\n# model %s \n#****** Result for:", model);        fprintf(ficlog,"\n# model %s \n#****** Result for:", model);
       for(j=1;j<=cptcoveff;j++){         /* It might not be a good idea to mix dummies and quantitative */
         printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);        /* for(j=1;j<=cptcoveff;j++){ /\* j=resultpos. Could be a loop on cptcovs: number of single dummy covariate in the result line as well as in the model *\/ */
         fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);        for(j=1;j<=cptcovs;j++){ /* j=resultpos. Could be a loop on cptcovs: number of single covariate (dummy or quantitative) in the result line as well as in the model */
         fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);          /* printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); /\* Output by variables in the resultline *\/ */
       }          /* Tvaraff[j] is the name of the dummy variable in position j in the equation model:
       for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */           * Tvaraff[1]@9={4, 3, 0, 0, 0, 0, 0, 0, 0}, in model=V5+V4+V3+V4*V3+V5*age
         printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);           * (V5 is quanti) V4 and V3 are dummies
         fprintf(ficrest," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);           * TnsdVar[4] is the position 1 and TnsdVar[3]=2 in codtabm(k,l)(V4  V3)=V4  V3
         fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);           *                                                              l=1 l=2
       }            *                                                           k=1  1   1   0   0
            *                                                           k=2  2   1   1   0
            *                                                           k=3 [1] [2]  0   1
            *                                                           k=4  2   2   1   1
            * If nres=1 result: V3=1 V4=0 then k=3 and outputs
            * If nres=2 result: V4=1 V3=0 then k=2 and outputs
            * nres=1 =>k=3 j=1 V4= nbcode[4][codtabm(3,1)=1)=0; j=2  V3= nbcode[3][codtabm(3,2)=2]=1
            * nres=2 =>k=2 j=1 V4= nbcode[4][codtabm(2,1)=2)=1; j=2  V3= nbcode[3][codtabm(2,2)=1]=0
            */
           /* 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  */
   /* 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(Dummy[modelresult[nres][j]]==0){/* Dummy variable of the variable in position modelresult in the model corresponding to j in resultline  */
             printf("V%d=%d ",Tvresult[nres][j],Tresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
             fprintf(ficlog,"V%d=%d ",Tvresult[nres][j],Tresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
             fprintf(ficrest,"V%d=%d ",Tvresult[nres][j],Tresult[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 */
               printf("fixed ");fprintf(ficlog,"fixed ");fprintf(ficrest,"fixed ");
             }else{
               printf("varyi ");fprintf(ficlog,"varyi ");fprintf(ficrest,"varyi ");
             }
             /* fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
             /* fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
           }else if(Dummy[modelresult[nres][j]]==1){ /* Quanti variable */
             /* For each selected (single) quantitative value */
             printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
             if(Fixed[modelresult[nres][j]]==0){ /* Fixed */
               printf("fixed ");fprintf(ficlog,"fixed ");fprintf(ficrest,"fixed ");
             }else{
               printf("varyi ");fprintf(ficlog,"varyi ");fprintf(ficrest,"varyi ");
             }
           }else{
             printf("Error in computing T_ Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=%d cptcovs=%d, cptcoveff=%d \n", nres, j, Dummy[modelresult[nres][j]],nres,j,modelresult[nres][j],cptcovs, cptcoveff);  /* end if dummy  or quanti */
             fprintf(ficlog,"Error in computing T_ Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=%d cptcovs=%d, cptcoveff=%d \n", nres, j, Dummy[modelresult[nres][j]],nres,j,modelresult[nres][j],cptcovs, cptcoveff);  /* end if dummy  or quanti */
             exit(1);
           }
         } /* End loop for each variable in the resultline */
         /* for (j=1; j<= nsq; j++){ /\* For each selected (single) quantitative value *\/ */
         /*        printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); /\* Wrong j is not in the equation model *\/ */
         /*        fprintf(ficrest," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */
         /*        fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */
         /* }       */
       fprintf(ficrest,"******\n");        fprintf(ficrest,"******\n");
       fprintf(ficlog,"******\n");        fprintf(ficlog,"******\n");
       printf("******\n");        printf("******\n");
Line 13713  Please run with mle=-1 to get a correct Line 13911  Please run with mle=-1 to get a correct
       fprintf(ficresstdeij,"\n#****** ");        fprintf(ficresstdeij,"\n#****** ");
       fprintf(ficrescveij,"\n#****** ");        fprintf(ficrescveij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++) {        for(j=1;j<=cptcoveff;j++) {
         fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);          fprintf(ficresstdeij,"V%d=%d ",Tvresult[nres][j],Tresult[nres][j]);
         fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);          /* fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
       }          /* fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
       for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */        }
         fprintf(ficresstdeij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);        for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value, TvarsQind gives the position of a quantitative in model equation  */
         fprintf(ficrescveij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);          fprintf(ficresstdeij," V%d=%f ",Tvar[TvarsQind[j]],Tqresult[nres][resultmodel[nres][TvarsQind[j]]]);
           fprintf(ficrescveij," V%d=%f ",Tvar[TvarsQind[j]],Tqresult[nres][resultmodel[nres][TvarsQind[j]]]);
       }         } 
       fprintf(ficresstdeij,"******\n");        fprintf(ficresstdeij,"******\n");
       fprintf(ficrescveij,"******\n");        fprintf(ficrescveij,"******\n");
Line 13726  Please run with mle=-1 to get a correct Line 13925  Please run with mle=-1 to get a correct
       fprintf(ficresvij,"\n#****** ");        fprintf(ficresvij,"\n#****** ");
       /* pstamp(ficresvij); */        /* pstamp(ficresvij); */
       for(j=1;j<=cptcoveff;j++)         for(j=1;j<=cptcoveff;j++) 
         fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[TnsdVar[Tvaraff[j]]])]);          fprintf(ficresvij,"V%d=%d ",Tvresult[nres][j],Tresult[nres][j]);
           /* fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[TnsdVar[Tvaraff[j]]])]); */
       for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */        for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
         /* fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); /\* To solve *\/ */          /* fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); /\* To solve *\/ */
         fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); /* Solved */          fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); /* Solved */
Line 13760  Please run with mle=-1 to get a correct Line 13960  Please run with mle=-1 to get a correct
           fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);            fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
         else          else
           fprintf(ficrest,"the age specific forward period (stable) prevalences in each health state \n");            fprintf(ficrest,"the age specific forward period (stable) prevalences in each health state \n");
         fprintf(ficrest,"# Age popbased mobilav e.. (std) ");          fprintf(ficrest,"# Age popbased mobilav e.. (std) "); /* Adding covariate values? */
         for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);          for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
         fprintf(ficrest,"\n");          fprintf(ficrest,"\n");
         /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */          /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
Line 13807  Please run with mle=-1 to get a correct Line 14007  Please run with mle=-1 to get a correct
       printf("done selection\n");fflush(stdout);        printf("done selection\n");fflush(stdout);
       fprintf(ficlog,"done selection\n");fflush(ficlog);        fprintf(ficlog,"done selection\n");fflush(ficlog);
               
     } /* End k selection */      } /* End k selection or end covariate selection for nres */
   
     printf("done State-specific expectancies\n");fflush(stdout);      printf("done State-specific expectancies\n");fflush(stdout);
     fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog);      fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog);
   
     /* variance-covariance of forward period prevalence*/      /* variance-covariance of forward period prevalence */
     varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff);      varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff);
   
           

Removed from v.1.333  
changed lines
  Added in v.1.336


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