* imach.c (Module): Many errors in graphs fixed with Vn*age covariates.

diff --git a/src/imach.c b/src/imach.c
index 39db19e53a038c605fcaa1f08347d91e263ffcd8..6710b39802ea88fb4cbd948e8d9fd39b9719b86c 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -1,6 +1,9 @@
 /* $Id$
   $State$
   $Log$
+  Revision 1.328  2022/07/27 17:40:48  brouard
+  Summary: valgrind bug fixed by initializing to zero DummyV as well as Tage
+
   Revision 1.327  2022/07/27 14:47:35  brouard
   Summary: Still a problem for one-step probabilities in case of quantitative variables
 
@@ -1189,7 +1192,7 @@ typedef struct {
 
 #define GNUPLOTPROGRAM "gnuplot"
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
-#define FILENAMELENGTH 132
+#define FILENAMELENGTH 256
 
 #define        GLOCK_ERROR_NOPATH              -1      /* empty path */
 #define        GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */
@@ -2821,8 +2824,8 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
     }
     for (k=1; k<=cptcovprod;k++){ /* For product without age */
       /* printf("prevalim Prod ij=%d k=%d  Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */
-      if(Dummy[Tvard[k][1]==0]){
-       if(Dummy[Tvard[k][2]==0]){
+      if(Dummy[Tvard[k][1]]==0){
+       if(Dummy[Tvard[k][2]]==0){
          cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
          /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
        }else{
@@ -2830,7 +2833,7 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
          /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; */
        }
       }else{
-       if(Dummy[Tvard[k][2]==0]){
+       if(Dummy[Tvard[k][2]]==0){
          cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
          /* cov[++k1]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; */
        }else{
@@ -3000,14 +3003,14 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
     }
     for (k=1; k<=cptcovprod;k++){ /* For product without age */
       /* printf("prevalim Prod ij=%d k=%d  Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */
-      if(Dummy[Tvard[k][1]==0]){
-       if(Dummy[Tvard[k][2]==0]){
+      if(Dummy[Tvard[k][1]]==0){
+       if(Dummy[Tvard[k][2]]==0){
          cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
        }else{
          cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k];
        }
       }else{
-       if(Dummy[Tvard[k][2]==0]){
+       if(Dummy[Tvard[k][2]]==0){
          cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
        }else{
          cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
@@ -3454,14 +3457,14 @@ double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, in
       for (k=1; k<=cptcovprod;k++){ /*  For product without age */
        /* printf("hPxij Prod ij=%d k=%d  Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */
        /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
-       if(Dummy[Tvard[k][1]==0]){
-         if(Dummy[Tvard[k][2]==0]){
+       if(Dummy[Tvard[k][1]]==0){
+         if(Dummy[Tvard[k][2]]==0){
            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
          }else{
            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k];
          }
        }else{
-         if(Dummy[Tvard[k][2]==0]){
+         if(Dummy[Tvard[k][2]]==0){
            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
          }else{
            cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
@@ -3575,14 +3578,14 @@ double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, do
       }
       for (k=1; k<=cptcovprod;k++){ /* Useless because included in cptcovn */
        cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
-       if(Dummy[Tvard[k][1]==0]){
-         if(Dummy[Tvard[k][2]==0]){
+       if(Dummy[Tvard[k][1]]==0){
+         if(Dummy[Tvard[k][2]]==0){
            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
          }else{
            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k];
          }
        }else{
-         if(Dummy[Tvard[k][2]==0]){
+         if(Dummy[Tvard[k][2]]==0){
            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
          }else{
            cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
@@ -5903,6 +5906,8 @@ void  concatwav(int wav[], int **dh, int **bh,  int **mw, int **s, double *agedc
 
 {
   /* Health expectancies, no variances */
+  /* cij is the combination in the list of combination of dummy covariates */
+  /* strstart is a string of time at start of computing */
   int i, j, nhstepm, hstepm, h, nstepm;
   int nhstepma, nstepma; /* Decreasing with age */
   double age, agelim, hf;
@@ -6979,8 +6984,8 @@ To be simple, these graphs help to understand the significativity of each parame
         /* 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]){
+        if(Dummy[Tvard[k][1]]==0){
+          if(Dummy[Tvard[k][2]]==0){
             cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,k)] * nbcode[Tvard[k][2]][codtabm(j1,k)];
             /* 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? */
@@ -6988,7 +6993,7 @@ To be simple, these graphs help to understand the significativity of each parame
             /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; */
           }
         }else{
-          if(Dummy[Tvard[k][2]==0]){
+          if(Dummy[Tvard[k][2]]==0){
             cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(j1,k)] * Tqinvresult[nres][Tvard[k][1]];
             /* cov[++k1]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; */
           }else{
@@ -7334,19 +7339,23 @@ divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3-%d.svg\">%s_%d-3-
 <img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); 
      /* Survival functions (period) in state j */
      for(cpt=1; cpt<=nlstate;cpt++){
-       fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
-<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
+       fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
+       fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
+       fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
      }
      /* State specific survival functions (period) */
      for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Survival functions in state %d and in any other live state (total).\
  And probability to be observed in various states (up to %d) being in state %d at different ages.      \
- <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> <img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
+ <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
+       fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
+       fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
      }
      /* Period (forward stable) prevalence in each health state */
      for(cpt=1; cpt<=nlstate;cpt++){
-       fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
-<img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
+       fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
+       fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"P_"),subdirf2(optionfilefiname,"P_"));
+      fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
      }
      if(prevbcast==1){
        /* Backward prevalence in each health state */
@@ -8407,41 +8416,49 @@ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar)
            /* for(j=3; j <=ncovmodel-nagesqr; j++) { */
            for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */
              /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
-             if(cptcovage >0){ /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
-               if(j==Tage[ij]) { /* Product by age  To be looked at!!*//* Bug valgrind */
-                 if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
-                   if(DummyV[j]==0){/* Bug valgrind */
-                     fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);;
-                   }else{ /* quantitative */
-                     fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */
-                     /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
+             switch(Typevar[j]){
+             case 1:
+               if(cptcovage >0){ /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
+                 if(j==Tage[ij]) { /* Product by age  To be looked at!!*//* Bug valgrind */
+                   if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
+                     if(DummyV[j]==0){/* Bug valgrind */
+                       fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);;
+                     }else{ /* quantitative */
+                       fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */
+                       /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
+                     }
+                     ij++;
                    }
-                   ij++;
                  }
-               } 
-             }else if(cptcovprod >0){
-               if(j==Tprod[ijp]) { /* */ 
-                 /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
-                 if(ijp <=cptcovprod) { /* Product */
-                   if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */
-                     if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */
-                       /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */
-                       fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]);
-                     }else{ /* Vn is dummy and Vm is quanti */
-                       /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
-                       fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
-                     }
-                   }else{ /* Vn*Vm Vn is quanti */
-                     if(DummyV[Tvard[ijp][2]]==0){
-                       fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]);
-                     }else{ /* Both quanti */
-                       fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+               }
+               break;
+             case 2:
+               if(cptcovprod >0){
+                 if(j==Tprod[ijp]) { /* */ 
+                   /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
+                   if(ijp <=cptcovprod) { /* Product */
+                     if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */
+                       if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */
+                         /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */
+                         fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]);
+                       }else{ /* Vn is dummy and Vm is quanti */
+                         /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
+                         fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+                       }
+                     }else{ /* Vn*Vm Vn is quanti */
+                       if(DummyV[Tvard[ijp][2]]==0){
+                         fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]);
+                       }else{ /* Both quanti */
+                         fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+                       }
                      }
+                     ijp++;
                    }
-                   ijp++;
-                 }
-               } /* end Tprod */
-             } else{  /* simple covariate */
+                 } /* end Tprod */
+               }
+               break;
+             case 0:
+               /* simple covariate */
                /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(k1,j)]); /\* Valgrind bug nbcode *\/ */
                if(Dummy[j]==0){
                  fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]); /*  */
@@ -8449,12 +8466,17 @@ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar)
                  fprintf(ficgp,"+p%d*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* */
                  /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
                }
-             } /* end simple */
+              /* end simple */
+               break;
+             default:
+               break;
+             } /* end switch */
            } /* end j */
-         }else{
-           i=i-ncovmodel;
-           if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
-             fprintf(ficgp," (1.");
+         }else{ /* k=k2 */
+           if(ng !=1 ){ /* For logit formula of log p11 is more difficult to get */
+             fprintf(ficgp," (1.");i=i-ncovmodel;
+           }else
+             i=i-ncovmodel;
          }
          
          if(ng != 1){
@@ -8467,17 +8489,78 @@ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar)
                fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(cpt-1)*ncovmodel,k3+(cpt-1)*ncovmodel+1,k3+(cpt-1)*ncovmodel+1+nagesqr);
               
              ij=1;
-             for(j=3; j <=ncovmodel-nagesqr; j++){
-                if(cptcovage >0){ 
-                  if((j-2)==Tage[ij]) { /* Bug valgrind */
-                    if(ij <=cptcovage) { /* Bug valgrind */
-                      fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);
-                      /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
-                      ij++;
-                    }
-                  }
-                }else
-                  fprintf(ficgp,"+p%d*%d",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);/* Valgrind bug nbcode */
+             ijp=1;
+             /* for(j=3; j <=ncovmodel-nagesqr; j++){ */
+             for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */
+               switch(Typevar[j]){
+               case 1:
+                 if(cptcovage >0){ 
+                   if(j==Tage[ij]) { /* Bug valgrind */
+                     if(ij <=cptcovage) { /* Bug valgrind */
+                       if(DummyV[j]==0){/* Bug valgrind */
+                         /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]); */
+                         /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,nbcode[Tvar[j]][codtabm(k1,j)]); */
+                         fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvar[j]]);
+                         /* fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);; */
+                         /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
+                       }else{ /* quantitative */
+                         /* fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /\* Tqinvresult in decoderesult *\/ */
+                         fprintf(ficgp,"+p%d*%f*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */
+                         /* fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /\* Tqinvresult in decoderesult *\/ */
+                         /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
+                       }
+                       ij++;
+                     }
+                   }
+                 }
+                 break;
+               case 2:
+                 if(cptcovprod >0){
+                   if(j==Tprod[ijp]) { /* */ 
+                     /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
+                     if(ijp <=cptcovprod) { /* Product */
+                       if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */
+                         if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */
+                           /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */
+                           fprintf(ficgp,"+p%d*%d*%d",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]);
+                           /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]); */
+                         }else{ /* Vn is dummy and Vm is quanti */
+                           /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
+                           fprintf(ficgp,"+p%d*%d*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+                           /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); */
+                         }
+                       }else{ /* Vn*Vm Vn is quanti */
+                         if(DummyV[Tvard[ijp][2]]==0){
+                           fprintf(ficgp,"+p%d*%d*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]);
+                           /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]); */
+                         }else{ /* Both quanti */
+                           fprintf(ficgp,"+p%d*%f*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+                           /* fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); */
+                         } 
+                       }
+                       ijp++;
+                     }
+                   } /* end Tprod */
+                 } /* end if */
+                 break;
+               case 0: 
+                 /* simple covariate */
+                 /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(k1,j)]); /\* Valgrind bug nbcode *\/ */
+                 if(Dummy[j]==0){
+                   /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]); /\*  *\/ */
+                   fprintf(ficgp,"+p%d*%d",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvar[j]]); /*  */
+                   /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]); /\*  *\/ */
+                 }else{ /* quantitative */
+                   fprintf(ficgp,"+p%d*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tqinvresult[nres][Tvar[j]]); /* */
+                   /* fprintf(ficgp,"+p%d*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /\* *\/ */
+                   /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
+                 }
+                 /* end simple */
+                 /* fprintf(ficgp,"+p%d*%d",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);/\* Valgrind bug nbcode *\/ */
+                 break;
+               default:
+                 break;
+               } /* end switch */
              }
              fprintf(ficgp,")");
            }
@@ -8486,7 +8569,7 @@ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar)
              fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"p%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
            else /* ng= 3 */
              fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"i%d%d\" ",  nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
-         }else{ /* end ng <> 1 */
+          }else{ /* end ng <> 1 */
            if( k !=k2) /* logit p11 is hard to draw */
              fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"logit(p%d%d)\" ",  nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
          }
@@ -12112,25 +12195,25 @@ Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numline
           * For k=4 covariates, h goes from 1 to m=2**k
           * codtabm(h,k)=  (1 & (h-1) >> (k-1)) + 1;
            * #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
-          *     h\k   1     2     3     4
-          *______________________________  
-          *     1 i=1 1 i=1 1 i=1 1 i=1 1
-          *     2     2     1     1     1
-          *     3 i=2 1     2     1     1
-          *     4     2     2     1     1
-          *     5 i=3 1 i=2 1     2     1
-          *     6     2     1     2     1
-          *     7 i=4 1     2     2     1
-          *     8     2     2     2     1
-          *     9 i=5 1 i=3 1 i=2 1     2
-          *    10     2     1     1     2
-          *    11 i=6 1     2     1     2
-          *    12     2     2     1     2
-          *    13 i=7 1 i=4 1     2     2    
-          *    14     2     1     2     2
-          *    15 i=8 1     2     2     2
-          *    16     2     2     2     2
-          */
+          *     h\k   1     2     3     4   *  h-1\k-1  4  3  2  1          
+          *______________________________   *______________________
+          *     1 i=1 1 i=1 1 i=1 1 i=1 1   *     0     0  0  0  0 
+          *     2     2     1     1     1   *     1     0  0  0  1 
+          *     3 i=2 1     2     1     1   *     2     0  0  1  0 
+          *     4     2     2     1     1   *     3     0  0  1  1 
+          *     5 i=3 1 i=2 1     2     1   *     4     0  1  0  0 
+          *     6     2     1     2     1   *     5     0  1  0  1 
+          *     7 i=4 1     2     2     1   *     6     0  1  1  0 
+          *     8     2     2     2     1   *     7     0  1  1  1 
+          *     9 i=5 1 i=3 1 i=2 1     2   *     8     1  0  0  0 
+          *    10     2     1     1     2   *     9     1  0  0  1 
+          *    11 i=6 1     2     1     2   *    10     1  0  1  0 
+          *    12     2     2     1     2   *    11     1  0  1  1 
+          *    13 i=7 1 i=4 1     2     2   *    12     1  1  0  0  
+          *    14     2     1     2     2   *    13     1  1  0  1 
+          *    15 i=8 1     2     2     2   *    14     1  1  1  0 
+          *    16     2     2     2     2   *    15     1  1  1  1          
+          */                                     
   /* How to do the opposite? From combination h (=1 to 2**k) how to get the value on the covariates? */
      /* from h=5 and m, we get then number of covariates k=log(m)/log(2)=4
      * and the value of each covariate?