imach/src/imach.c - diff

Return to imach.c CVS log

Up to [Local Repository] / imach / src

Diff for /imach/src/imach.c between versions 1.236 and 1.240

-version 1.236, 2016/08/25 10:50:18
+version 1.240, 2016/08/29 07:53:18
  Line 1
  /* $Id$
    $State$
    $Log$
+   Revision 1.240  2016/08/29 07:53:18  brouard
+   Summary: Better
+   Revision 1.239  2016/08/26 15:51:03  brouard
+   Summary: Improvement in Powell output in order to copy and paste
+   Author:
+   Revision 1.238  2016/08/26 14:23:35  brouard
+   Summary: Starting tests of 0.99
+   Revision 1.237  2016/08/26 09:20:19  brouard
+   Summary: to valgrind
    Revision 1.236  2016/08/25 10:50:18  brouard
    *** empty log message ***
- Line 1026  double dval;
+ Line 1040  double dval;
  #define FTOL 1.0e-10
  #define NRANSI
  #define ITMAX 200
+ #define ITPOWMAX 20 /* This is now multiplied by the number of parameters */
  #define TOL 2.0e-4
- Line 1114  int *TvarsQind;
+ Line 1129  int *TvarsQind;
  #define MAXRESULTLINES 10
  int nresult=0;
  int TKresult[MAXRESULTLINES];
- double Tresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
+ int Tresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
+ int Tinvresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
  int Tvresult[MAXRESULTLINES][NCOVMAX]; /* For dummy variable , variable # (output) */
  double Tqresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , value (output) */
+ double Tqinvresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , value (output) */
  int Tvqresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , variable # (output) */
  /* int *TDvar; /\**< TDvar[1]=4,  TDvarF[2]=3, TDvar[3]=6  in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 *\/ */
- Line 1140  double *Tvalsel; /**< Selected modality
+ Line 1157  double *Tvalsel; /**< Selected modality
  int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product */
  int *Fixed; /** Fixed[k] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */
  int *Dummy; /** Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */
+ int *DummyV; /** Dummy[v] 0=dummy (0 1), 1 quantitative */
+ int *FixedV; /** FixedV[v] 0 fixed, 1 varying */
  int *Tage;
  int anyvaryingduminmodel=0; /**< Any varying dummy in Model=1 yes, 0 no, to avoid a loop on waves in freq */
  int *Tmodelind; /** Tmodelind[Tvaraff[3]]=9 for V1 position,Tvaraff[1]@9={4, 3, 1, 0, 0, 0, 0, 0, 0}, model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/
- Line 1149  int *Ndum; /** Freq of modality (tricode
+ Line 1168  int *Ndum; /** Freq of modality (tricode
  /* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */
  int **Tvard;
  int *Tprod;/**< Gives the k position of the k1 product */
+ /* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3  */
  int *Tposprod; /**< Gives the k1 product from the k position */
- /* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3
+    /* if  V2+V1+V1*V4+age*V3+V3*V2   TProd[k1=2]=5 (V3*V2) */
-    if  V2+V1+V1*V4+age*V3+V3*V2   TProd[k1=2]=5 (V3*V2)
+    /* Tposprod[k]=k1 , Tposprod[3]=1, Tposprod[5(V3*V2)]=2 (2nd product without age) */
-    Tposprod[k]=k1 , Tposprod[3]=1, Tposprod[5]=2
- */
  int cptcovprod, *Tvaraff, *invalidvarcomb;
  double *lsurv, *lpop, *tpop;
- Line 2055  void powell(double p[], double **xi, int
+ Line 2073  void powell(double p[], double **xi, int
   void linmin(double p[], double xi[], int n, double *fret,
                                                   double (*func)(double []),int *flat);
  #endif
-   int i,ibig,j;
+  int i,ibig,j,jk,k;
    double del,t,*pt,*ptt,*xit;
    double directest;
    double fp,fptt;
- Line 2087  void powell(double p[], double **xi, int
+ Line 2105  void powell(double p[], double **xi, int
      fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
  /*     fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */
      for (i=1;i<=n;i++) {
-       printf(" %d %.12f",i, p[i]);
-       fprintf(ficlog," %d %.12lf",i, p[i]);
        fprintf(ficrespow," %.12lf", p[i]);
      }
+     fprintf(ficrespow,"\n");fflush(ficrespow);
+     printf("\n#model=  1      +     age ");
+     fprintf(ficlog,"\n#model=  1      +     age ");
+     if(nagesqr==1){
+         printf("  + age*age  ",Tvar[j]);
+         fprintf(ficlog,"  + age*age  ",Tvar[j]);
+     }
+     for(j=1;j <=ncovmodel-2;j++){
+       if(Typevar[j]==0) {
+         printf("  +      V%d  ",Tvar[j]);
+         fprintf(ficlog,"  +      V%d  ",Tvar[j]);
+       }else if(Typevar[j]==1) {
+         printf("  +    V%d*age ",Tvar[j]);
+         fprintf(ficlog,"  +    V%d*age ",Tvar[j]);
+       }else if(Typevar[j]==2) {
+         printf("  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+         fprintf(ficlog,"  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+       }
+     }
      printf("\n");
+ /*     printf("12   47.0114589    0.0154322   33.2424412    0.3279905    2.3731903  */
+ /* 13  -21.5392400    0.1118147    1.2680506    1.2973408   -1.0663662  */
      fprintf(ficlog,"\n");
-     fprintf(ficrespow,"\n");fflush(ficrespow);
+     for(i=1,jk=1; i <=nlstate; i++){
+       for(k=1; k <=(nlstate+ndeath); k++){
+         if (k != i) {
+           printf("%d%d ",i,k);
+           fprintf(ficlog,"%d%d ",i,k);
+           for(j=1; j <=ncovmodel; j++){
+             printf("%12.7f ",p[jk]);
+             fprintf(ficlog,"%12.7f ",p[jk]);
+             jk++;
+           }
+           printf("\n");
+           fprintf(ficlog,"\n");
+         }
+       }
+     }
      if(*iter <=3){
        tml = *localtime(&rcurr_time);
        strcpy(strcurr,asctime(&tml));
- Line 2211  void powell(double p[], double **xi, int
+ Line 2262  void powell(double p[], double **xi, int
        free_vector(pt,1,n);
        return;
      } /* enough precision */
-     if (*iter == ITMAX) nrerror("powell exceeding maximum iterations.");
+     if (*iter == ITMAX*n) nrerror("powell exceeding maximum iterations.");
      for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */
        ptt[j]=2.0*p[j]-pt[j];
        xit[j]=p[j]-pt[j];
- Line 2410  void powell(double p[], double **xi, int
+ Line 2461  void powell(double p[], double **xi, int
        cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];
        /* printf("prevalim Quantitative k=%d  TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
      }
-     for (k=1; k<=cptcovage;k++){
+     for (k=1; k<=cptcovage;k++){  /* For product with age */
        if(Dummy[Tvar[Tage[k]]]){
          cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
        } else{
- Line 2418  void powell(double p[], double **xi, int
+ Line 2469  void powell(double p[], double **xi, int
        }
        /* printf("prevalim Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
      }
-     for (k=1; k<=cptcovprod;k++){ /*  */
+     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]); */
-       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]){
+           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]){
+           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]];
+         }
+       }
      }
      /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
      /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
- Line 4073  void  freqsummary(char fileres[], int ia
+ Line 4136  void  freqsummary(char fileres[], int ia
      fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
      exit(0);
    }
    strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));
    if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {
      printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
- Line 4083  void  freqsummary(char fileres[], int ia
+ Line 4146  void  freqsummary(char fileres[], int ia
    }
    else{
      fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%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",\
              fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
    }
-   fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\n",fileresphtm, fileresphtm);
+   fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm);
    strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
    if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
      printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
      fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
      fflush(ficlog);
      exit(70);
-   }
+   } else{
-   else{
      fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%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",\
              fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
    }
-   fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
+   fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
    freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
    j1=0;
    /* j=ncoveff;  /\* Only fixed dummy covariates *\/ */
    j=cptcoveff;  /* Only dummy covariates of the model */
    if (cptcovn<1) {j=1;ncodemax[1]=1;}
    first=1;
    /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels:
       reference=low_education V1=0,V2=0
       med_educ                V1=1 V2=0,
       high_educ               V1=0 V2=1
       Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff
    */
    for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives V4=0, V3=0 for example, fixed or varying covariates */
      posproptt=0.;
      /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
        scanf("%d", i);*/
      for (i=-5; i<=nlstate+ndeath; i++)
        for (jk=-5; jk<=nlstate+ndeath; jk++)
-                                 for(m=iagemin; m <= iagemax+3; m++)
+         for(m=iagemin; m <= iagemax+3; m++)
-                                         freq[i][jk][m]=0;
+           freq[i][jk][m]=0;
      for (i=1; i<=nlstate; i++)  {
        for(m=iagemin; m <= iagemax+3; m++)
-                                 prop[i][m]=0;
+         prop[i][m]=0;
        posprop[i]=0;
        pospropt[i]=0;
      }
- Line 4141  Title=%s <br>Datafile=%s Firstpass=%d La
+ Line 4203  Title=%s <br>Datafile=%s Firstpass=%d La
      /*  meanqt[m][z1]=0.; */
      /*   } */
      /* } */
      dateintsum=0;
      k2cpt=0;
      /* For that combination of covariate j1, we count and print the frequencies in one pass */
- Line 4220  Title=%s <br>Datafile=%s Firstpass=%d La
+ Line 4282  Title=%s <br>Datafile=%s Firstpass=%d La
      } /* end iind = 1 to imx */
      /* prop[s][age] is feeded for any initial and valid live state as well as
         freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */
      /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
      pstamp(ficresp);
-     /* if  (ncoveff>0) { */
+     if  (cptcoveff>0){
-     if  (cptcoveff>0) {
        fprintf(ficresp, "\n#********** Variable ");
        fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");
        fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");
+       fprintf(ficlog, "\n#********** Variable ");
        for (z1=1; z1<=cptcoveff; z1++){
-         fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+         if(DummyV[z1]){
-         fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+           fprintf(ficresp, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-         fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+           fprintf(ficresphtm, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+           fprintf(ficresphtmfr, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+           fprintf(ficlog, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+         }else{
+           fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+           fprintf(ficresphtm, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+           fprintf(ficresphtmfr, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+           fprintf(ficlog, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+         }
        }
        fprintf(ficresp, "**********\n#");
        fprintf(ficresphtm, "**********</h3>\n");
        fprintf(ficresphtmfr, "**********</h3>\n");
-       fprintf(ficlog, "\n#********** Variable ");
-       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
        fprintf(ficlog, "**********\n");
      }
      fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
      for(i=1; i<=nlstate;i++) {
-       fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
+       fprintf(ficresp, " Age 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(ficresp, "\n");
      fprintf(ficresphtm, "\n");
      /* Header of frequency table by age */
      fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");
      fprintf(ficresphtmfr,"<th>Age</th> ");
- Line 4259  Title=%s <br>Datafile=%s Firstpass=%d La
+ Line 4327  Title=%s <br>Datafile=%s Firstpass=%d La
        }
      }
      fprintf(ficresphtmfr, "\n");
      /* For each age */
      for(iage=iagemin; iage <= iagemax+3; iage++){
        fprintf(ficresphtm,"<tr>");
        if(iage==iagemax+1){
-                                 fprintf(ficlog,"1");
+         fprintf(ficlog,"1");
-                                 fprintf(ficresphtmfr,"<tr><th>0</th> ");
+         fprintf(ficresphtmfr,"<tr><th>0</th> ");
        }else if(iage==iagemax+2){
-                                 fprintf(ficlog,"0");
+         fprintf(ficlog,"0");
-                                 fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
+         fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
        }else if(iage==iagemax+3){
-                                 fprintf(ficlog,"Total");
+         fprintf(ficlog,"Total");
-                                 fprintf(ficresphtmfr,"<tr><th>Total</th> ");
+         fprintf(ficresphtmfr,"<tr><th>Total</th> ");
        }else{
-                                 if(first==1){
+         if(first==1){
-                                         first=0;
+           first=0;
-                                         printf("See log file for details...\n");
+           printf("See log file for details...\n");
-                                 }
+         }
-                                 fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage);
+         fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage);
-                                 fprintf(ficlog,"Age %d", iage);
+         fprintf(ficlog,"Age %d", iage);
        }
        for(jk=1; jk <=nlstate ; jk++){
-                                 for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
+         for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
-                                         pp[jk] += freq[jk][m][iage];
+           pp[jk] += freq[jk][m][iage];
        }
        for(jk=1; jk <=nlstate ; jk++){
-                                 for(m=-1, pos=0; m <=0 ; m++)
+         for(m=-1, pos=0; m <=0 ; m++)
-                                         pos += freq[jk][m][iage];
+           pos += freq[jk][m][iage];
-                                 if(pp[jk]>=1.e-10){
+         if(pp[jk]>=1.e-10){
-                                         if(first==1){
+           if(first==1){
-                                                 printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
+             printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
-                                         }
+           }
-                                         fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
+           fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
-                                 }else{
+         }else{
-                                         if(first==1)
+           if(first==1)
-                                                 printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
+             printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
-                                         fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
+           fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
-                                 }
+         }
        }
        for(jk=1; jk <=nlstate ; jk++){
-                                 /* posprop[jk]=0; */
+         /* posprop[jk]=0; */
-                                 for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */
+         for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */
-                                         pp[jk] += freq[jk][m][iage];
+           pp[jk] += freq[jk][m][iage];
        } /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */
        for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){
-                                 pos += pp[jk]; /* pos is the total number of transitions until this age */
+         pos += pp[jk]; /* pos is the total number of transitions until this age */
-                                 posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state
+         posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state
-                                                                                                                                                                         from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
+                                           from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
-                                 pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state
+         pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state
-                                                                                                                                                                 from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
+                                         from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
        }
        for(jk=1; jk <=nlstate ; jk++){
-                                 if(pos>=1.e-5){
+         if(pos>=1.e-5){
-                                         if(first==1)
+           if(first==1)
-                                                 printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
+             printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
-                                         fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
+           fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
-                                 }else{
+         }else{
-                                         if(first==1)
+           if(first==1)
-                                                 printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
+             printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
-                                         fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
+           fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
-                                 }
+         }
-                                 if( iage <= iagemax){
+         if( iage <= iagemax){
-                                         if(pos>=1.e-5){
+           if(pos>=1.e-5){
-                                                 fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
+             fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
-                                                 fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
+             fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
-                                                 /*probs[iage][jk][j1]= pp[jk]/pos;*/
+             /*probs[iage][jk][j1]= pp[jk]/pos;*/
-                                                 /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/
+             /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/
-                                         }
+           }
-                                         else{
+           else{
-                                                 fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);
+             fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);
-                                                 fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);
+             fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);
-                                         }
+           }
-                                 }
+         }
-                                 pospropt[jk] +=posprop[jk];
+         pospropt[jk] +=posprop[jk];
        } /* end loop jk */
        /* pospropt=0.; */
        for(jk=-1; jk <=nlstate+ndeath; jk++){
-                                 for(m=-1; m <=nlstate+ndeath; m++){
+         for(m=-1; m <=nlstate+ndeath; m++){
-                                         if(freq[jk][m][iage] !=0 ) { /* minimizing output */
+           if(freq[jk][m][iage] !=0 ) { /* minimizing output */
-                                                 if(first==1){
+             if(first==1){
-                                                         printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
+               printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
-                                                 }
+             }
-                                                 fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
+             fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
-                                         }
+           }
-                                         if(jk!=0 && m!=0)
+           if(jk!=0 && m!=0)
-                                                 fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]);
+             fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]);
-                                 }
+         }
        } /* end loop jk */
        posproptt=0.;
        for(jk=1; jk <=nlstate; jk++){
-                                 posproptt += pospropt[jk];
+         posproptt += pospropt[jk];
        }
        fprintf(ficresphtmfr,"</tr>\n ");
        if(iage <= iagemax){
-                                 fprintf(ficresp,"\n");
+         fprintf(ficresp,"\n");
-                                 fprintf(ficresphtm,"</tr>\n");
+         fprintf(ficresphtm,"</tr>\n");
        }
        if(first==1)
-                                 printf("Others in log...\n");
+         printf("Others in log...\n");
        fprintf(ficlog,"\n");
      } /* end loop age iage */
      fprintf(ficresphtm,"<tr><th>Tot</th>");
      for(jk=1; jk <=nlstate ; jk++){
        if(posproptt < 1.e-5){
-                                 fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt);
+         fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt);
        }else{
-                                 fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);
+         fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);
        }
      }
      fprintf(ficresphtm,"</tr>\n");
- Line 4385  Title=%s <br>Datafile=%s Firstpass=%d La
+ Line 4453  Title=%s <br>Datafile=%s Firstpass=%d La
      fprintf(ficresphtmfr,"</table>\n");
    } /* end selected combination of covariate j1 */
    dateintmean=dateintsum/k2cpt;
    fclose(ficresp);
    fclose(ficresphtm);
    fclose(ficresphtmfr);
- Line 4901  void  concatwav(int wav[], int **dh, int
+ Line 4969  void  concatwav(int wav[], int **dh, int
    double ***p3mat;
    double eip;
-   pstamp(ficreseij);
+   /* pstamp(ficreseij); */
    fprintf(ficreseij,"# (a) Life expectancies by health status at initial age and (b) health expectancies by health status at initial age\n");
    fprintf(ficreseij,"# Age");
    for(i=1; i<=nlstate;i++){
- Line 5266  void  concatwav(int wav[], int **dh, int
+ Line 5334  void  concatwav(int wav[], int **dh, int
     fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
     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,"# Selected quantitative variables and dummies");
+    for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+      fprintf(ficresprobmorprev," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+    }
+    for(j=1;j<=cptcoveff;j++)
+      fprintf(ficresprobmorprev,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(ij,j)]);
+    fprintf(ficresprobmorprev,"\n");
     fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
     for(j=nlstate+1; j<=(nlstate+ndeath);j++){
       fprintf(ficresprobmorprev," p.%-d SE",j);
- Line 5986  void printinghtml(char fileresu[], char
+ Line 6062  void printinghtml(char fileresu[], char
                    int popforecast, int prevfcast, int backcast, int estepm , \
                    double jprev1, double mprev1,double anprev1, double dateprev1, \
                    double jprev2, double mprev2,double anprev2, double dateprev2){
-   int jj1, k1, i1, cpt;
+   int jj1, k1, i1, cpt, k4, nres;
     fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
     <li><a href='#secondorder'>Result files (second order (variance)</a>\n \
  </ul>");
+    fprintf(fichtm,"<ul><li> model=1+age+%s\n \
+ </ul>", model);
     fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n");
     fprintf(fichtm,"<li>- Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file)<br/>\n",
             jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm"));
- Line 6025  void printinghtml(char fileresu[], char
+ Line 6103  void printinghtml(char fileresu[], char
     if (cptcovn < 1) {m=1;ncodemax[1]=1;}
     jj1=0;
+    for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k1=1; k1<=m;k1++){
+      if(TKresult[nres]!= k1)
+        continue;
       /* for(i1=1; i1<=ncodemax[k1];i1++){ */
       jj1++;
       if (cptcovn > 0) {
         fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
         for (cpt=1; cpt<=cptcoveff;cpt++){
-          fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
+          fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],(int)Tresult[nres][cpt]);
-          printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout);
+          printf(" V%d=%d ",Tvresult[nres][cpt],Tresult[nres][cpt]);fflush(stdout);
+          /* fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); */
+          /* printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); */
         }
+        for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+         fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+         printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);fflush(stdout);
+       }
         /* if(nqfveff+nqtveff 0) */ /* Test to be done */
         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
         if(invalidvarcomb[k1]){
- Line 6146  See page 'Matrix of variance-covariance
+ Line 6235  See page 'Matrix of variance-covariance
     if (cptcovn < 1) {m=1;ncodemax[1]=1;}
     jj1=0;
+    for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k1=1; k1<=m;k1++){
+      if(TKresult[nres]!= k1)
+        continue;
       /* for(i1=1; i1<=ncodemax[k1];i1++){ */
       jj1++;
       if (cptcovn > 0) {
         fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
         for (cpt=1; cpt<=cptcoveff;cpt++)  /**< cptcoveff number of variables */
-          fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
+          fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],Tresult[nres][cpt]);
+          /* fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); */
+        for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+         fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+       }
         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
         if(invalidvarcomb[k1]){
- Line 6182  void printinggnuplot(char fileresu[], ch
+ Line 6280  void printinggnuplot(char fileresu[], ch
    char dirfileres[132],optfileres[132];
    char gplotcondition[132];
-   int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,ij=0,l=0;
+   int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,ij=0, ijp=0, l=0;
    int lv=0, vlv=0, kl=0;
    int ng=0;
    int vpopbased;
- Line 6231  void printinggnuplot(char fileresu[], ch
+ Line 6329  void printinggnuplot(char fileresu[], ch
    strcpy(dirfileres,optionfilefiname);
    strcpy(optfileres,"vpl");
    /* 1eme*/
-   for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */
+   for (cpt=1; cpt<= nlstate ; cpt ++){ /* For each live state */
-     for (k1=1; k1<= m ; k1 ++) /* For each valid combination of covariate */
+     for (k1=1; k1<= m ; k1 ++){ /* For each valid combination of covariate */
        for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-     /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
+         /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
-       if(TKresult[nres]!= k1)
+         if(TKresult[nres]!= k1)
-         continue;
+           continue;
-       /* We are interested in selected combination by the resultline */
+         /* We are interested in selected combination by the resultline */
-       printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt);
+         printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt);
-       fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files  and live state =%d ", cpt);
+         fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files  and live state =%d ", cpt);
-       for (k=1; k<=cptcoveff; k++){    /* For each covariate k get corresponding value lv for combination k1 */
+         for (k=1; k<=cptcoveff; k++){    /* For each covariate k get corresponding value lv for combination k1 */
-         lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
+           lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
-         /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+           /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-         /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+           /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-         /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+           /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-         vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */
+           vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */
-         /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
+           /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
-         printf(" V%d=%d ",Tvaraff[k],vlv);
+           printf(" V%d=%d ",Tvaraff[k],vlv);
-         fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+           fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-       }
+         }
-       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+         for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-         printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+           printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-         fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+           fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-       }
+         }
-       printf("\n#\n");
+         printf("\n#\n");
-       fprintf(ficgp,"\n#\n");
+         fprintf(ficgp,"\n#\n");
-       if(invalidvarcomb[k1]){
+         if(invalidvarcomb[k1]){
-         fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+           fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
-         continue;
+           continue;
-       }
+         }
-       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
+         fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
-       fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
+         fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
-       fprintf(ficgp,"set xlabel \"Age\" \n\
+         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
- set ylabel \"Probability\" \n             \
- set ter svg size 640, 480\n                                             \
+         for (i=1; i<= nlstate ; i ++) {
- plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+           if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+           else        fprintf(ficgp," %%*lf (%%*lf)");
-       for (i=1; i<= nlstate ; i ++) {
+         }
-         if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+         fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
-         else        fprintf(ficgp," %%*lf (%%*lf)");
+         for (i=1; i<= nlstate ; i ++) {
-       }
+           if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
-       fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+           else fprintf(ficgp," %%*lf (%%*lf)");
-       for (i=1; i<= nlstate ; i ++) {
+         }
-         if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+         fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
-         else fprintf(ficgp," %%*lf (%%*lf)");
+         for (i=1; i<= nlstate ; i ++) {
-       }
+           if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
-       fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+           else fprintf(ficgp," %%*lf (%%*lf)");
-       for (i=1; i<= nlstate ; i ++) {
+         }
-         if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+         fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
-         else fprintf(ficgp," %%*lf (%%*lf)");
+         if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */
-       }
+           /* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */
-       fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
+           fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */
-       if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */
+           if(cptcoveff ==0){
-         /* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */
+             fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ",        2+(cpt-1),  cpt );
-         fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */
+           }else{
-         if(cptcoveff ==0){
+             kl=0;
-           fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ",  2+(cpt-1),  cpt );
+             for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */
-         }else{
+               lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
-           kl=0;
+               /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-           for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */
+               /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-             lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+               /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-             /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+               vlv= nbcode[Tvaraff[k]][lv];
-             /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-             /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-             vlv= nbcode[Tvaraff[k]][lv];
-             kl++;
-             /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
-             /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */
-             /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */
-             /* ''  u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
-             if(k==cptcoveff){
-               fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
-+(cpt-1),  cpt );  /* 4 or 6 ?*/
-             }else{
-               fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
                kl++;
-             }
+               /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
-           } /* end covariate */
+               /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */
-         } /* end if no covariate */
+               /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */
-       } /* end if backcast */
+               /* ''  u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
-       fprintf(ficgp,"\nset out \n");
+               if(k==cptcoveff){
+                 fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
++(cpt-1),  cpt );  /* 4 or 6 ?*/
+               }else{
+                 fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
+                 kl++;
+               }
+             } /* end covariate */
+           } /* end if no covariate */
+         } /* end if backcast */
+         fprintf(ficgp,"\nset out \n");
+       } /* nres */
      } /* k1 */
    } /* cpt */
    /*2 eme*/
-   for (k1=1; k1<= m ; k1 ++)
+   for (k1=1; k1<= m ; k1 ++){
-   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-     if(TKresult[nres]!= k1)
+       if(TKresult[nres]!= k1)
-       continue;
+         continue;
-     fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
+       fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
-     for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-       lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
-       /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-       /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-       /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-       vlv= nbcode[Tvaraff[k]][lv];
-       fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-     }
-     for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-       printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-       fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-     }
-     fprintf(ficgp,"\n#\n");
-     if(invalidvarcomb[k1]){
-       fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
-       continue;
-     }
-     fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
-     for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
-       if(vpopbased==0)
-         fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
-       else
-         fprintf(ficgp,"\nreplot ");
-       for (i=1; i<= nlstate+1 ; i ++) {
-         k=2*i;
-         fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
-         for (j=1; j<= nlstate+1 ; j ++) {
-           if (j==i) fprintf(ficgp," %%lf (%%lf)");
-           else fprintf(ficgp," %%*lf (%%*lf)");
-         }
-         if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
-         else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
-         fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
-         for (j=1; j<= nlstate+1 ; j ++) {
-           if (j==i) fprintf(ficgp," %%lf (%%lf)");
-           else fprintf(ficgp," %%*lf (%%*lf)");
-         }
-         fprintf(ficgp,"\" t\"\" w l lt 0,");
-         fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
-         for (j=1; j<= nlstate+1 ; j ++) {
-           if (j==i) fprintf(ficgp," %%lf (%%lf)");
-           else fprintf(ficgp," %%*lf (%%*lf)");
-         }
-         if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
-         else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
-       } /* state */
-     } /* vpopbased */
-     fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
-   } /* k1 end 2 eme*/
-   /*3eme*/
-   for (k1=1; k1<= m ; k1 ++)
-   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-     if(TKresult[nres]!= k)
-       continue;
-     for (cpt=1; cpt<= nlstate ; cpt ++) {
-       fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files:  combination=%d state=%d",k1, cpt);
        for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
          lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
          /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
- Line 6390  plot [%.f:%.f] \"%s\" every :::%d::%d u
+ Line 6426  plot [%.f:%.f] \"%s\" every :::%d::%d u
          vlv= nbcode[Tvaraff[k]][lv];
          fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
        }
+       /* for(k=1; k <= ncovds; k++){ */
        for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+         printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
          fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
        }
        fprintf(ficgp,"\n#\n");
        if(invalidvarcomb[k1]){
          fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
          continue;
        }
-       /*       k=2+nlstate*(2*cpt-2); */
+       fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
-       k=2+(nlstate+1)*(cpt-1);
+       for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
-       fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
+         if(vpopbased==0)
-       fprintf(ficgp,"set ter svg size 640, 480\n\
+           fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
+         else
+           fprintf(ficgp,"\nreplot ");
+         for (i=1; i<= nlstate+1 ; i ++) {
+           k=2*i;
+           fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
+           for (j=1; j<= nlstate+1 ; j ++) {
+             if (j==i) fprintf(ficgp," %%lf (%%lf)");
+             else fprintf(ficgp," %%*lf (%%*lf)");
+           }
+           if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
+           else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
+           fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+           for (j=1; j<= nlstate+1 ; j ++) {
+             if (j==i) fprintf(ficgp," %%lf (%%lf)");
+             else fprintf(ficgp," %%*lf (%%*lf)");
+           }
+           fprintf(ficgp,"\" t\"\" w l lt 0,");
+           fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+           for (j=1; j<= nlstate+1 ; j ++) {
+             if (j==i) fprintf(ficgp," %%lf (%%lf)");
+             else fprintf(ficgp," %%*lf (%%*lf)");
+           }
+           if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
+           else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
+         } /* state */
+       } /* vpopbased */
+       fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
+     } /* end nres */
+   } /* k1 end 2 eme*/
+   /*3eme*/
+   for (k1=1; k1<= m ; k1 ++){
+     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+       if(TKresult[nres]!= k1)
+         continue;
+       for (cpt=1; cpt<= nlstate ; cpt ++) {
+         fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files:  combination=%d state=%d",k1, cpt);
+         for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+           lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+           /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+           /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+           /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+           vlv= nbcode[Tvaraff[k]][lv];
+           fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+         }
+         for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+           fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+         }
+         fprintf(ficgp,"\n#\n");
+         if(invalidvarcomb[k1]){
+           fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+           continue;
+         }
+         /*       k=2+nlstate*(2*cpt-2); */
+         k=2+(nlstate+1)*(cpt-1);
+         fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
+         fprintf(ficgp,"set ter svg size 640, 480\n\
  plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt);
-       /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
+         /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
-         for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+           for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
-         fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+           fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
-         fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
+           fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
-         for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+           for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
-         fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+           fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
-       */
+         */
-       for (i=1; i< nlstate ; i ++) {
+         for (i=1; i< nlstate ; i ++) {
-         fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
+           fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
-         /*      fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
+           /*    fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
-       }
+         }
-       fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
+         fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
-     }
+       }
-   }
+     } /* end nres */
+   } /* end kl 3eme */
    /* 4eme */
    /* Survival functions (period) from state i in state j by initial state i */
-   for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+   for (k1=1; k1<=m; k1++){    /* For each covariate and each value */
-   for(nres=1; nres <= nresult; nres++) /* For each resultline */
+     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-     if(TKresult[nres]!= k)
+       if(TKresult[nres]!= k1)
-       continue;
-     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
-       fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
-         lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
-         /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-         /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-         /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-         vlv= nbcode[Tvaraff[k]][lv];
-         fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-       }
-       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-         fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-       }
-       fprintf(ficgp,"\n#\n");
-       if(invalidvarcomb[k1]){
-         fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
          continue;
-       }
+       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state cpt*/
+         fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
-       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
+         for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
+           lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
- set ter svg size 640, 480\n                                                                                                                                                                                     \
+           /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
- unset log y\n                                                                                                                                                                                                                                           \
+           /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
- plot [%.f:%.f]  ", ageminpar, agemaxpar);
+           /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-       k=3;
+           vlv= nbcode[Tvaraff[k]][lv];
-       for (i=1; i<= nlstate ; i ++){
+           fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-         if(i==1){
-           fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
-         }else{
-           fprintf(ficgp,", '' ");
          }
-         l=(nlstate+ndeath)*(i-1)+1;
+         for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-         fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
+           fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-         for (j=2; j<= nlstate+ndeath ; j ++)
+         }
-           fprintf(ficgp,"+$%d",k+l+j-1);
+         fprintf(ficgp,"\n#\n");
-         fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt);
+         if(invalidvarcomb[k1]){
-       } /* nlstate */
+           fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
-       fprintf(ficgp,"\nset out\n");
+           continue;
-     } /* end cpt state*/
+         }
-   } /* end covariate */
+         fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
+         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
+ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
+         k=3;
+         for (i=1; i<= nlstate ; i ++){
+           if(i==1){
+             fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+           }else{
+             fprintf(ficgp,", '' ");
+           }
+           l=(nlstate+ndeath)*(i-1)+1;
+           fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
+           for (j=2; j<= nlstate+ndeath ; j ++)
+             fprintf(ficgp,"+$%d",k+l+j-1);
+           fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt);
+         } /* nlstate */
+         fprintf(ficgp,"\nset out\n");
+       } /* end cpt state*/
+     } /* end nres */
+   } /* end covariate k1 */
  /* 5eme */
    /* Survival functions (period) from state i in state j by final state j */
-   for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
+   for (k1=1; k1<= m ; k1++){ /* For each covariate combination if any */
-   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-     if(TKresult[nres]!= k1)
+       if(TKresult[nres]!= k1)
-       continue;
-     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
-       fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
-       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-         lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
-         /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
-         /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
-         /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-         vlv= nbcode[Tvaraff[k]][lv];
-         fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
-       }
-       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-         fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-       }
-       fprintf(ficgp,"\n#\n");
-       if(invalidvarcomb[k1]){
-         fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
          continue;
-       }
+       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
+         fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
+         for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+           lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+           /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+           /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+           /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+           vlv= nbcode[Tvaraff[k]][lv];
+           fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+         }
+         for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+           fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+         }
+         fprintf(ficgp,"\n#\n");
+         if(invalidvarcomb[k1]){
+           fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+           continue;
+         }
-       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
+         fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
-       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
+         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
- set ter svg size 640, 480\n                                             \
+ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
- unset log y\n                                                           \
+         k=3;
- plot [%.f:%.f]  ", ageminpar, agemaxpar);
+         for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
-       k=3;
+           if(j==1)
-       for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
+             fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
-         if(j==1)
+           else
-           fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+             fprintf(ficgp,", '' ");
-         else
+           l=(nlstate+ndeath)*(cpt-1) +j;
-           fprintf(ficgp,", '' ");
+           fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l);
-         l=(nlstate+ndeath)*(cpt-1) +j;
+           /* for (i=2; i<= nlstate+ndeath ; i ++) */
-         fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l);
+           /*   fprintf(ficgp,"+$%d",k+l+i-1); */
-         /* for (i=2; i<= nlstate+ndeath ; i ++) */
+           fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j);
-         /*   fprintf(ficgp,"+$%d",k+l+i-1); */
+         } /* nlstate */
-         fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j);
+         fprintf(ficgp,", '' ");
-       } /* nlstate */
+         fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1);
-       fprintf(ficgp,", '' ");
+         for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
-       fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1);
+           l=(nlstate+ndeath)*(cpt-1) +j;
-       for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
+           if(j < nlstate)
-         l=(nlstate+ndeath)*(cpt-1) +j;
+             fprintf(ficgp,"$%d +",k+l);
-         if(j < nlstate)
+           else
-           fprintf(ficgp,"$%d +",k+l);
+             fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);
-         else
+         }
-           fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);
+         fprintf(ficgp,"\nset out\n");
-       }
+       } /* end cpt state*/
-       fprintf(ficgp,"\nset out\n");
+     } /* end covariate */
-     } /* end cpt state*/
+   } /* end nres */
-   } /* end covariate */
  /* 6eme */
    /* CV preval stable (period) for each covariate */
-   for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
+   for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
+   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+     if(TKresult[nres]!= k1)
+       continue;
      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
        fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
- Line 6538  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6637  plot [%.f:%.f]  ", ageminpar, agemaxpar)
          vlv= nbcode[Tvaraff[k]][lv];
          fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
        }
+       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+         fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+       }
        fprintf(ficgp,"\n#\n");
        if(invalidvarcomb[k1]){
          fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- Line 6546  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6648  plot [%.f:%.f]  ", ageminpar, agemaxpar)
        fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);
        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
- set ter svg size 640, 480\n                                             \
+ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
- unset log y\n                                                           \
- plot [%.f:%.f]  ", ageminpar, agemaxpar);
        k=3; /* Offset */
        for (i=1; i<= nlstate ; i ++){
          if(i==1)
- Line 6569  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6669  plot [%.f:%.f]  ", ageminpar, agemaxpar)
  /* 7eme */
    if(backcast == 1){
      /* CV back preval stable (period) for each covariate */
-     for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
+     for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
+     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+       if(TKresult[nres]!= k1)
+         continue;
        for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
          fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
          for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
- Line 6580  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6683  plot [%.f:%.f]  ", ageminpar, agemaxpar)
            vlv= nbcode[Tvaraff[k]][lv];
            fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
          }
+         for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+           fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+         }
          fprintf(ficgp,"\n#\n");
          if(invalidvarcomb[k1]){
            fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- Line 6588  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6694  plot [%.f:%.f]  ", ageminpar, agemaxpar)
          fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
          fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
- set ter svg size 640, 480\n                                             \
+ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
- unset log y\n                                                           \
- plot [%.f:%.f]  ", ageminpar, agemaxpar);
          k=3; /* Offset */
          for (i=1; i<= nlstate ; i ++){
            if(i==1)
- Line 6616  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6720  plot [%.f:%.f]  ", ageminpar, agemaxpar)
    if(prevfcast==1){
      /* Projection from cross-sectional to stable (period) for each covariate */
-     for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
+     for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
+     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+       if(TKresult[nres]!= k1)
+         continue;
        for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
          fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
          for (k=1; k<=cptcoveff; k++){    /* For each correspondig covariate value  */
- Line 6627  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6734  plot [%.f:%.f]  ", ageminpar, agemaxpar)
            vlv= nbcode[Tvaraff[k]][lv];
            fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
          }
+         for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+           fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+         }
          fprintf(ficgp,"\n#\n");
          if(invalidvarcomb[k1]){
            fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
- Line 6636  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6746  plot [%.f:%.f]  ", ageminpar, agemaxpar)
          fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
          fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
          fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
- set ter svg size 640, 480\n                                             \
+ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
- unset log y\n                                                           \
- plot [%.f:%.f]  ", ageminpar, agemaxpar);
          for (i=1; i<= nlstate+1 ; i ++){  /* nlstate +1 p11 p21 p.1 */
            /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
            /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */
- Line 6705  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6813  plot [%.f:%.f]  ", ageminpar, agemaxpar)
    } /* End if prevfcast */
-   /* proba elementaires */
+   /* 9eme writing MLE parameters */
-   fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
+   fprintf(ficgp,"\n##############\n#9eme MLE estimated parameters\n#############\n");
    for(i=1,jk=1; i <=nlstate; i++){
      fprintf(ficgp,"# initial state %d\n",i);
      for(k=1; k <=(nlstate+ndeath); k++){
- Line 6723  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6831  plot [%.f:%.f]  ", ageminpar, agemaxpar)
    fprintf(ficgp,"##############\n#\n");
    /*goto avoid;*/
-   fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n");
+   /* 10eme Graphics of probabilities or incidences using written MLE parameters */
+   fprintf(ficgp,"\n##############\n#10eme Graphics of probabilities or incidences\n#############\n");
    fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n");
    fprintf(ficgp,"# logi(p12/p11)=p1 +p2*age +p3*age*age+ p4*V1+ p5*V1*age\n");
    fprintf(ficgp,"# logi(p13/p11)=a13+b13*age+c13age*age+d13*V1+e13*V1*age\n");
- Line 6738  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6847  plot [%.f:%.f]  ", ageminpar, agemaxpar)
    fprintf(ficgp,"#       +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n");
    fprintf(ficgp,"#\n");
    for(ng=1; ng<=3;ng++){ /* Number of graphics: first is logit, 2nd is probabilities, third is incidences per year*/
-     fprintf(ficgp,"# ng=%d\n",ng);
+     fprintf(ficgp,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year\n");
-     fprintf(ficgp,"#   jk=1 to 2^%d=%d\n",cptcoveff,m);
+     fprintf(ficgp,"#model=%s \n",model);
-     for(jk=1; jk <=m; jk++) {
+     fprintf(ficgp,"# Type of graphic ng=%d\n",ng);
-       fprintf(ficgp,"#    jk=%d\n",jk);
+     fprintf(ficgp,"#   jk=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */
+     for(jk=1; jk <=m; jk++)  /* For each combination of covariate */
+     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+       if(TKresult[nres]!= jk)
+         continue;
+       fprintf(ficgp,"# Combination of dummy  jk=%d and ",jk);
+       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+         fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+       }
+       fprintf(ficgp,"\n#\n");
        fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng);
        fprintf(ficgp,"\nset ter svg size 640, 480 ");
        if (ng==1){
- Line 6782  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6900  plot [%.f:%.f]  ", ageminpar, agemaxpar)
                break;
              }
              ij=1;/* To be checked else nbcode[0][0] wrong */
-             for(j=3; j <=ncovmodel-nagesqr; j++) {
+             ijp=1; /* product no age */
+             /* 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(ij <=cptcovage) { /* Bug valgrind */
+               if(j==Tage[ij]) { /* Product by age */
-                 if((j-2)==Tage[ij]) { /* Bug valgrind */
+                 if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
-                   fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+                   if(DummyV[j]==0){
-                   /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+                     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(jk,Tvar[j-2])]); */
+                   }
                    ij++;
                  }
-               }
+               }else if(j==Tprod[ijp]) { /* */
-               else
+                 /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
-                 fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); /* Valgrind bug nbcode */
+                 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(jk,j)],nbcode[Tvard[ijp][2]][codtabm(jk,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(jk,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++;
+                 }
+               } else{  /* simple covariate */
+                 /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(jk,j)]); /\* Valgrind bug nbcode *\/ */
+                 if(Dummy[j]==0){
+                   fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]); /*  */
+                 }else{ /* quantitative */
+                   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(jk,Tvar[j-2])]); */
+                 }
+               } /* end simple */
+             } /* end j */
            }else{
              i=i-ncovmodel;
              if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
- Line 6811  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6961  plot [%.f:%.f]  ", ageminpar, agemaxpar)
                ij=1;
                for(j=3; j <=ncovmodel-nagesqr; j++){
-                 if(ij <=cptcovage) { /* Bug valgrind */
+                 if((j-2)==Tage[ij]) { /* Bug valgrind */
-                   if((j-2)==Tage[ij]) { /* Bug valgrind */
+                   if(ij <=cptcovage) { /* Bug valgrind */
                      fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
                      /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
                      ij++;
- Line 7076  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 7226  plot [%.f:%.f]  ", ageminpar, agemaxpar)
        fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
      }
      for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-       printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+       fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-       fprintf(ficlog," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
      }
      fprintf(ficresf," yearproj age");
      for(j=1; j<=nlstate+ndeath;j++){
- Line 7686  int readdata(char datafile[], int firsto
+ Line 7835  int readdata(char datafile[], int firsto
    /*-------- data file ----------*/
    FILE *fic;
    char dummy[]="                         ";
-   int i=0, j=0, n=0, iv=0;
+   int i=0, j=0, n=0, iv=0, v;
    int lstra;
    int linei, month, year,iout;
    char line[MAXLINE], linetmp[MAXLINE];
    char stra[MAXLINE], strb[MAXLINE];
    char *stratrunc;
+   DummyV=ivector(1,NCOVMAX); /* 1 to 3 */
+   FixedV=ivector(1,NCOVMAX); /* 1 to 3 */
+   for(v=1; v <=ncovcol;v++){
+     DummyV[v]=0;
+     FixedV[v]=0;
+   }
+   for(v=ncovcol+1; v <=ncovcol+nqv;v++){
+     DummyV[v]=1;
+     FixedV[v]=0;
+   }
+   for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){
+     DummyV[v]=0;
+     FixedV[v]=1;
+   }
+   for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){
+     DummyV[v]=1;
+     FixedV[v]=1;
+   }
+   for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){
+     printf("Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
+     fprintf(ficlog,"Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
+   }
    if((fic=fopen(datafile,"r"))==NULL)    {
      printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);
- Line 7723  int readdata(char datafile[], int firsto
+ Line 7894  int readdata(char datafile[], int firsto
      /* Loops on waves */
      for (j=maxwav;j>=1;j--){
        for (iv=nqtv;iv>=1;iv--){  /* Loop  on time varying quantitative variables */
-                                 cutv(stra, strb, line, ' ');
+         cutv(stra, strb, line, ' ');
-                                 if(strb[0]=='.') { /* Missing value */
+         if(strb[0]=='.') { /* Missing value */
-                                         lval=-1;
+           lval=-1;
-                                         cotqvar[j][iv][i]=-1; /* 0.0/0.0 */
+           cotqvar[j][iv][i]=-1; /* 0.0/0.0 */
-                                         cotvar[j][ntv+iv][i]=-1; /* For performance reasons */
+           cotvar[j][ntv+iv][i]=-1; /* For performance reasons */
-                                         if(isalpha(strb[1])) { /* .m or .d Really Missing value */
+           if(isalpha(strb[1])) { /* .m or .d Really Missing value */
-                                                 printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);
+             printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);
-                                                 fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);fflush(ficlog);
+             fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);fflush(ficlog);
-                                                 return 1;
+             return 1;
-                                         }
+           }
-                                 }else{
+         }else{
-                                         errno=0;
+           errno=0;
-                                         /* what_kind_of_number(strb); */
+           /* what_kind_of_number(strb); */
-                                         dval=strtod(strb,&endptr);
+           dval=strtod(strb,&endptr);
-                                         /* if( strb[0]=='\0' || (*endptr != '\0')){ */
+           /* if( strb[0]=='\0' || (*endptr != '\0')){ */
-                                         /* if(strb != endptr && *endptr == '\0') */
+           /* if(strb != endptr && *endptr == '\0') */
-                                         /*    dval=dlval; */
+           /*    dval=dlval; */
-                                         /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */
+           /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */
-                                         if( strb[0]=='\0' || (*endptr != '\0')){
+           if( strb[0]=='\0' || (*endptr != '\0')){
-                                                 printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, nqtv, j,maxwav);
+             printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, nqtv, j,maxwav);
-                                                 fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line, iv, nqtv, j,maxwav);fflush(ficlog);
+             fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line, iv, nqtv, j,maxwav);fflush(ficlog);
-                                                 return 1;
+             return 1;
-                                         }
+           }
-                                         cotqvar[j][iv][i]=dval;
+           cotqvar[j][iv][i]=dval;
-                                         cotvar[j][ntv+iv][i]=dval;
+           cotvar[j][ntv+iv][i]=dval;
-                                 }
+         }
-                                 strcpy(line,stra);
+         strcpy(line,stra);
        }/* end loop ntqv */
        for (iv=ntv;iv>=1;iv--){  /* Loop  on time varying dummies */
-                                 cutv(stra, strb, line, ' ');
+         cutv(stra, strb, line, ' ');
-                                 if(strb[0]=='.') { /* Missing value */
+         if(strb[0]=='.') { /* Missing value */
-                                         lval=-1;
+           lval=-1;
-                                 }else{
+         }else{
-                                         errno=0;
+           errno=0;
-                                         lval=strtol(strb,&endptr,10);
+           lval=strtol(strb,&endptr,10);
-                                         /*      if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
+           /*    if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
-                                         if( strb[0]=='\0' || (*endptr != '\0')){
+           if( strb[0]=='\0' || (*endptr != '\0')){
-                                                 printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th dummy covariate out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, ntv, j,maxwav);
+             printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th dummy covariate out of %d measured at wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, ntv, j,maxwav);
-                                                 fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d dummy covariate out of %d measured wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, ntv,j,maxwav);fflush(ficlog);
+             fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d dummy covariate out of %d measured wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,iv, ntv,j,maxwav);fflush(ficlog);
-                                                 return 1;
+             return 1;
-                                         }
+           }
-                                 }
+         }
-                                 if(lval <-1 || lval >1){
+         if(lval <-1 || lval >1){
-                                         printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
+           printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
   Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
   for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
-  For example, for multinomial values like 1, 2 and 3,\n                                                                 \
+  For example, for multinomial values like 1, 2 and 3,\n                 \
-  build V1=0 V2=0 for the reference value (1),\n                                                                                                 \
+  build V1=0 V2=0 for the reference value (1),\n                         \
-         V1=1 V2=0 for (2) \n                                                                                                                                                                            \
+         V1=1 V2=0 for (2) \n                                            \
   and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
-  output of IMaCh is often meaningless.\n                                                                                                                                \
+  output of IMaCh is often meaningless.\n                                \
   Exiting.\n",lval,linei, i,line,j);
-                                         fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
+           fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
   Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
   for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
-  For example, for multinomial values like 1, 2 and 3,\n                                                                 \
+  For example, for multinomial values like 1, 2 and 3,\n                 \
-  build V1=0 V2=0 for the reference value (1),\n                                                                                                 \
+  build V1=0 V2=0 for the reference value (1),\n                         \
-         V1=1 V2=0 for (2) \n                                                                                                                                                                            \
+         V1=1 V2=0 for (2) \n                                            \
   and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
-  output of IMaCh is often meaningless.\n                                                                                                                                \
+  output of IMaCh is often meaningless.\n                                \
   Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
-                                         return 1;
+           return 1;
-                                 }
+         }
-                                 cotvar[j][iv][i]=(double)(lval);
+         cotvar[j][iv][i]=(double)(lval);
-                                 strcpy(line,stra);
+         strcpy(line,stra);
        }/* end loop ntv */
        /* Statuses  at wave */
        cutv(stra, strb, line, ' ');
        if(strb[0]=='.') { /* Missing value */
-                                 lval=-1;
+         lval=-1;
        }else{
-                                 errno=0;
+         errno=0;
-                                 lval=strtol(strb,&endptr,10);
+         lval=strtol(strb,&endptr,10);
-                                 /*      if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
+         /*      if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
-                                 if( strb[0]=='\0' || (*endptr != '\0')){
+         if( strb[0]=='\0' || (*endptr != '\0')){
-                                         printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);
+           printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);
-                                         fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
+           fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
-                                         return 1;
+           return 1;
-                                 }
+         }
        }
        s[j][i]=lval;
- Line 8020  int decoderesult ( char resultline[], in
+ Line 8191  int decoderesult ( char resultline[], in
      if(Typevar[k1]==0){ /* Single covariate in model */
        match=0;
        for(k2=1; k2 <=j;k2++){/* result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
-         if(Tvar[k1]==Tvarsel[k2]) {/* Tvar[2]=5 == Tvarsel[1]=4   */
+         if(Tvar[k1]==Tvarsel[k2]) {/* Tvar[1]=5 == Tvarsel[2]=5   */
            modelresult[k2]=k1;/* modelresult[2]=1 modelresult[1]=2  modelresult[3]=3  modelresult[6]=4 modelresult[9]=5 */
            match=1;
            break;
- Line 8036  int decoderesult ( char resultline[], in
+ Line 8207  int decoderesult ( char resultline[], in
      match=0;
      for(k1=1; k1<= cptcovt ;k1++){ /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
        if(Typevar[k1]==0){ /* Single */
-         if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=5 == Tvarsel[1]=4   */
+         if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=4 == Tvarsel[1]=4   */
            resultmodel[k1]=k2;  /* resultmodel[2]=1 resultmodel[1]=2  resultmodel[3]=3  resultmodel[6]=4 resultmodel[9]=5 */
            ++match;
          }
- Line 8063  int decoderesult ( char resultline[], in
+ Line 8234  int decoderesult ( char resultline[], in
    /*    6 1 0 1 */ /* V4=1, V3=0, V1=1 */
    /*    7 0 1 1 */
    /*    8 1 1 1 */
+   /* V(Tvresult)=Tresult V4=1 V3=0 V1=1 Tresult[nres=1][2]=0 */
+   /* V(Tvqresult)=Tqresult V5=25.1 V2=8 Tqresult[nres=1][1]=25.1 */
+   /* V5*age V5 known which value for nres?  */
+   /* Tqinvresult[2]=8 Tqinvresult[1]=25.1  */
    for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* model line */
      if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Single dummy */
-       k3= resultmodel[k1]; /* resultmodel[2] = 1=k3 */
+       k3= resultmodel[k1]; /* resultmodel[2(V4)] = 1=k3 */
        k2=(int)Tvarsel[k3]; /*  Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */
        k+=Tvalsel[k3]*pow(2,k4);  /*  Tvalsel[1]=1  */
-       Tresult[nres][k4+1]=Tvalsel[k3];
+       Tresult[nres][k4+1]=Tvalsel[k3];/* Tresult[nres][1]=1(V4=1)  Tresult[nres][2]=0(V3=0) */
-       Tvresult[nres][k4+1]=(int)Tvarsel[k3];
+       Tvresult[nres][k4+1]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */
+       Tinvresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */
        printf("Decoderesult Dummy k=%d, V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k2, k3, (int)Tvalsel[k3], k4);
        k4++;;
      }  else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Single quantitative */
        k3q= resultmodel[k1]; /* resultmodel[2] = 1=k3 */
        k2q=(int)Tvarsel[k3q]; /*  Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */
-       Tqresult[nres][k4q+1]=Tvalsel[k3q];
+       Tqresult[nres][k4q+1]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */
-       Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q];
+       Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */
+       Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */
        printf("Decoderesult Quantitative nres=%d, V(k2q=V%d)= Tvalsel[%d]=%d, Tvarsel[%d]=%f\n",nres, k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]);
        k4q++;;
      }
- Line 8102  int decodemodel( char model[], int lasto
+ Line 8279  int decodemodel( char model[], int lasto
          * - Tvard[k]  p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 .
          */
  {
-   int i, j, k, ks;
+   int i, j, k, ks, v;
    int  j1, k1, k2, k3, k4;
    char modelsav[80];
    char stra[80], strb[80], strc[80], strd[80],stre[80];
- Line 8308  Dummy[k] 0=dummy (0 1), 1 quantitative (
+ Line 8485  Dummy[k] 0=dummy (0 1), 1 quantitative (
  Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product \n\
  Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
  Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
+   for(k=1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;}
    for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */
      if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
        Fixed[k]= 0;
- Line 8333  Dummy[k] 0=dummy (0 1), 1 quantitative (
+ Line 8510  Dummy[k] 0=dummy (0 1), 1 quantitative (
        TvarFind[ncovf]=k;
        TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
        TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
-     }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/ /* Only simple fixed quantitative variable */
+     }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){/* Remind that product Vn*Vm are added in k Only simple fixed quantitative variable */
        Fixed[k]= 0;
        Dummy[k]= 1;
        nqfveff++;
- Line 8386  Dummy[k] 0=dummy (0 1), 1 quantitative (
+ Line 8563  Dummy[k] 0=dummy (0 1), 1 quantitative (
        TvarA[ncova]=Tvar[k];
        TvarAind[ncova]=k;
        if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */
          Fixed[k]= 2;
          Dummy[k]= 2;
          modell[k].maintype= ATYPE;
          modell[k].subtype= APFD;
          /* ncoveff++; */
        }else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/
          Fixed[k]= 2;
          Dummy[k]= 3;
          modell[k].maintype= ATYPE;
          modell[k].subtype= APFQ;                /*      Product age * fixed quantitative */
          /* nqfveff++;  /\* Only simple fixed quantitative variable *\/ */
        }else if( Tvar[k] <=ncovcol+nqv+ntv ){
          Fixed[k]= 3;
          Dummy[k]= 2;
          modell[k].maintype= ATYPE;
          modell[k].subtype= APVD;                /*      Product age * varying dummy */
          /* ntveff++; /\* Only simple time varying dummy variable *\/ */
        }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){
          Fixed[k]= 3;
          Dummy[k]= 3;
          modell[k].maintype= ATYPE;
          modell[k].subtype= APVQ;                /*      Product age * varying quantitative */
          /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */
        }
      }else if (Typevar[k] == 2) {  /* product without age */
        k1=Tposprod[k];
        if(Tvard[k1][1] <=ncovcol){
          if(Tvard[k1][2] <=ncovcol){
            Fixed[k]= 1;
            Dummy[k]= 0;
            modell[k].maintype= FTYPE;
            modell[k].subtype= FPDD;              /*      Product fixed dummy * fixed dummy */
            ncovf++; /* Fixed variables without age */
            TvarF[ncovf]=Tvar[k];
            TvarFind[ncovf]=k;
          }else if(Tvard[k1][2] <=ncovcol+nqv){
            Fixed[k]= 0;  /* or 2 ?*/
            Dummy[k]= 1;
            modell[k].maintype= FTYPE;
            modell[k].subtype= FPDQ;              /*      Product fixed dummy * fixed quantitative */
            ncovf++; /* Varying variables without age */
            TvarF[ncovf]=Tvar[k];
            TvarFind[ncovf]=k;
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
            Fixed[k]= 1;
            Dummy[k]= 0;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDD;              /*      Product fixed dummy * varying dummy */
            ncovv++; /* Varying variables without age */
            TvarV[ncovv]=Tvar[k];
            TvarVind[ncovv]=k;
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
            Fixed[k]= 1;
            Dummy[k]= 1;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDQ;              /*      Product fixed dummy * varying quantitative */
            ncovv++; /* Varying variables without age */
            TvarV[ncovv]=Tvar[k];
            TvarVind[ncovv]=k;
          }
        }else if(Tvard[k1][1] <=ncovcol+nqv){
          if(Tvard[k1][2] <=ncovcol){
            Fixed[k]= 0;  /* or 2 ?*/
            Dummy[k]= 1;
            modell[k].maintype= FTYPE;
            modell[k].subtype= FPDQ;              /*      Product fixed quantitative * fixed dummy */
            ncovf++; /* Fixed variables without age */
            TvarF[ncovf]=Tvar[k];
            TvarFind[ncovf]=k;
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
            Fixed[k]= 1;
            Dummy[k]= 1;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDQ;              /*      Product fixed quantitative * varying dummy */
            ncovv++; /* Varying variables without age */
            TvarV[ncovv]=Tvar[k];
            TvarVind[ncovv]=k;
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
            Fixed[k]= 1;
            Dummy[k]= 1;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPQQ;              /*      Product fixed quantitative * varying quantitative */
            ncovv++; /* Varying variables without age */
            TvarV[ncovv]=Tvar[k];
            TvarVind[ncovv]=k;
            ncovv++; /* Varying variables without age */
            TvarV[ncovv]=Tvar[k];
            TvarVind[ncovv]=k;
          }
        }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){
          if(Tvard[k1][2] <=ncovcol){
            Fixed[k]= 1;
            Dummy[k]= 1;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDD;              /*      Product time varying dummy * fixed dummy */
            ncovv++; /* Varying variables without age */
            TvarV[ncovv]=Tvar[k];
            TvarVind[ncovv]=k;
          }else if(Tvard[k1][2] <=ncovcol+nqv){
            Fixed[k]= 1;
            Dummy[k]= 1;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDQ;              /*      Product time varying dummy * fixed quantitative */
            ncovv++; /* Varying variables without age */
            TvarV[ncovv]=Tvar[k];
            TvarVind[ncovv]=k;
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
            Fixed[k]= 1;
            Dummy[k]= 0;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDD;              /*      Product time varying dummy * time varying dummy */
            ncovv++; /* Varying variables without age */
            TvarV[ncovv]=Tvar[k];
            TvarVind[ncovv]=k;
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
            Fixed[k]= 1;
            Dummy[k]= 1;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDQ;              /*      Product time varying dummy * time varying quantitative */
            ncovv++; /* Varying variables without age */
            TvarV[ncovv]=Tvar[k];
            TvarVind[ncovv]=k;
          }
        }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){
          if(Tvard[k1][2] <=ncovcol){
            Fixed[k]= 1;
            Dummy[k]= 1;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDQ;              /*      Product time varying quantitative * fixed dummy */
            ncovv++; /* Varying variables without age */
            TvarV[ncovv]=Tvar[k];
            TvarVind[ncovv]=k;
          }else if(Tvard[k1][2] <=ncovcol+nqv){
            Fixed[k]= 1;
            Dummy[k]= 1;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPQQ;              /*      Product time varying quantitative * fixed quantitative */
            ncovv++; /* Varying variables without age */
            TvarV[ncovv]=Tvar[k];
            TvarVind[ncovv]=k;
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
            Fixed[k]= 1;
            Dummy[k]= 1;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDQ;              /*      Product time varying quantitative * time varying dummy */
            ncovv++; /* Varying variables without age */
            TvarV[ncovv]=Tvar[k];
            TvarVind[ncovv]=k;
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
            Fixed[k]= 1;
            Dummy[k]= 1;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPQQ;              /*      Product time varying quantitative * time varying quantitative */
            ncovv++; /* Varying variables without age */
            TvarV[ncovv]=Tvar[k];
            TvarVind[ncovv]=k;
          }
        }else{
          printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
          fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
-       } /* end k1 */
+       } /*end k1*/
      }else{
        printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
        fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
- Line 8923  int prevalence_limit(double *p, double *
+ Line 9100  int prevalence_limit(double *p, double *
    i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
    if (cptcovn < 1){i1=1;}
-   for(nres=1; nres <= nresult; nres++) /* For each resultline */
+   for(k=1; k<=i1;k++){ /* For each combination k of dummy covariates in the model */
-   for(k=1; k<=i1;k++){
+     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-     if(TKresult[nres]!= k)
+       if(TKresult[nres]!= k)
-       continue;
+         continue;
-   /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
+       /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
-     /* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
+       /* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
-     //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+       //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
-     /* k=k+1; */
+       /* k=k+1; */
-     /* to clean */
+       /* to clean */
-     //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+       //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
-     fprintf(ficrespl,"#******");
+       fprintf(ficrespl,"#******");
-     printf("#******");
+       printf("#******");
-     fprintf(ficlog,"#******");
+       fprintf(ficlog,"#******");
-     for(j=1;j<=cptcoveff ;j++) {/* all covariates */
+       for(j=1;j<=cptcoveff ;j++) {/* all covariates */
-       fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); /* Here problem for varying dummy*/
+         fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); /* Here problem for varying dummy*/
-       printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+         printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-       fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+         fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-     }
+       }
-     for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-       printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+         printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-       fprintf(ficlog," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+         fprintf(ficrespl," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-     }
+         fprintf(ficlog," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
-     fprintf(ficrespl,"******\n");
+       }
-     printf("******\n");
+       fprintf(ficrespl,"******\n");
-     fprintf(ficlog,"******\n");
+       printf("******\n");
-     if(invalidvarcomb[k]){
+       fprintf(ficlog,"******\n");
-       printf("\nCombination (%d) ignored because no case \n",k);
+       if(invalidvarcomb[k]){
-       fprintf(ficrespl,"#Combination (%d) ignored because no case \n",k);
+         printf("\nCombination (%d) ignored because no case \n",k);
-       fprintf(ficlog,"\nCombination (%d) ignored because no case \n",k);
+         fprintf(ficrespl,"#Combination (%d) ignored because no case \n",k);
-                                                 continue;
+         fprintf(ficlog,"\nCombination (%d) ignored because no case \n",k);
-     }
+         continue;
+       }
-     fprintf(ficrespl,"#Age ");
+       fprintf(ficrespl,"#Age ");
-     for(j=1;j<=cptcoveff;j++) {
+       for(j=1;j<=cptcoveff;j++) {
-       fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+         fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-     }
+       }
-     for(i=1; i<=nlstate;i++) fprintf(ficrespl,"  %d-%d   ",i,i);
+       for(i=1; i<=nlstate;i++) fprintf(ficrespl,"  %d-%d   ",i,i);
-     fprintf(ficrespl,"Total Years_to_converge\n");
+       fprintf(ficrespl,"Total Years_to_converge\n");
-     for (age=agebase; age<=agelim; age++){
+       for (age=agebase; age<=agelim; age++){
-       /* for (age=agebase; age<=agebase; age++){ */
+         /* for (age=agebase; age<=agebase; age++){ */
-       prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k, nres);
+         prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k, nres);
-       fprintf(ficrespl,"%.0f ",age );
+         fprintf(ficrespl,"%.0f ",age );
-       for(j=1;j<=cptcoveff;j++)
+         for(j=1;j<=cptcoveff;j++)
-         fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+           fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-       tot=0.;
+         tot=0.;
-       for(i=1; i<=nlstate;i++){
+         for(i=1; i<=nlstate;i++){
-         tot +=  prlim[i][i];
+           tot +=  prlim[i][i];
-         fprintf(ficrespl," %.5f", prlim[i][i]);
+           fprintf(ficrespl," %.5f", prlim[i][i]);
-       }
+         }
-       fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
+         fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
-     } /* Age */
+       } /* Age */
-     /* was end of cptcod */
+       /* was end of cptcod */
-   } /* cptcov */
+     } /* cptcov */
+   } /* nres */
    return 0;
  }
- Line 9018  int back_prevalence_limit(double *p, dou
+ Line 9197  int back_prevalence_limit(double *p, dou
    i1=pow(2,cptcoveff);
    if (cptcovn < 1){i1=1;}
-   for(nres=1; nres <= nresult; nres++) /* For each resultline */
+   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-   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 */
-     if(TKresult[nres]!= k)
+       if(TKresult[nres]!= k)
-       continue;
+         continue;
-     //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+       //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
-     fprintf(ficresplb,"#******");
+       fprintf(ficresplb,"#******");
-     printf("#******");
+       printf("#******");
-     fprintf(ficlog,"#******");
+       fprintf(ficlog,"#******");
-     for(j=1;j<=cptcoveff ;j++) {/* all covariates */
+       for(j=1;j<=cptcoveff ;j++) {/* all covariates */
-       fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+         fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-       printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+         printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-       fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+         fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,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][j]);
+         printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-       fprintf(ficresplb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+         fprintf(ficresplb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-       fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+         fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
-     }
+       }
-     fprintf(ficresplb,"******\n");
+       fprintf(ficresplb,"******\n");
-     printf("******\n");
+       printf("******\n");
-     fprintf(ficlog,"******\n");
+       fprintf(ficlog,"******\n");
-     if(invalidvarcomb[k]){
+       if(invalidvarcomb[k]){
-       printf("\nCombination (%d) ignored because no cases \n",k);
+         printf("\nCombination (%d) ignored because no cases \n",k);
-       fprintf(ficresplb,"#Combination (%d) ignored because no cases \n",k);
+         fprintf(ficresplb,"#Combination (%d) ignored because no cases \n",k);
-       fprintf(ficlog,"\nCombination (%d) ignored because no cases \n",k);
+         fprintf(ficlog,"\nCombination (%d) ignored because no cases \n",k);
-       continue;
+         continue;
-     }
+       }
-     fprintf(ficresplb,"#Age ");
+       fprintf(ficresplb,"#Age ");
-     for(j=1;j<=cptcoveff;j++) {
+       for(j=1;j<=cptcoveff;j++) {
-       fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+         fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-     }
+       }
-     for(i=1; i<=nlstate;i++) fprintf(ficresplb,"  %d-%d   ",i,i);
+       for(i=1; i<=nlstate;i++) fprintf(ficresplb,"  %d-%d   ",i,i);
-     fprintf(ficresplb,"Total Years_to_converge\n");
+       fprintf(ficresplb,"Total Years_to_converge\n");
-     for (age=agebase; age<=agelim; age++){
+       for (age=agebase; age<=agelim; age++){
-       /* for (age=agebase; age<=agebase; age++){ */
+         /* for (age=agebase; age<=agebase; age++){ */
-       if(mobilavproj > 0){
+         if(mobilavproj > 0){
-         /* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */
+           /* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */
-         /* bprevalim(bprlim, mobaverage, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
+           /* bprevalim(bprlim, mobaverage, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
-         bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k);
+           bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k);
-       }else if (mobilavproj == 0){
+         }else if (mobilavproj == 0){
-         printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+           printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
-         fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+           fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
-         exit(1);
+           exit(1);
-       }else{
+         }else{
-         /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
+           /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
-         bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);
+           bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);
-       }
+         }
-       fprintf(ficresplb,"%.0f ",age );
+         fprintf(ficresplb,"%.0f ",age );
-       for(j=1;j<=cptcoveff;j++)
+         for(j=1;j<=cptcoveff;j++)
-         fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+           fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-       tot=0.;
+         tot=0.;
-       for(i=1; i<=nlstate;i++){
+         for(i=1; i<=nlstate;i++){
-         tot +=  bprlim[i][i];
+           tot +=  bprlim[i][i];
-         fprintf(ficresplb," %.5f", bprlim[i][i]);
+           fprintf(ficresplb," %.5f", bprlim[i][i]);
-       }
+         }
-       fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
+         fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
-     } /* Age */
+       } /* Age */
-     /* was end of cptcod */
+       /* was end of cptcod */
-   } /* cptcov */
+     } /* end of any combination */
+   } /* end of nres */
    /* hBijx(p, bage, fage); */
    /* fclose(ficrespijb); */
- Line 9172  int hPijx(double *p, int bage, int fage)
+ Line 9351  int hPijx(double *p, int bage, int fage)
          int ageminl;
    int hstepm;
    int nhstepm;
-   int h, i, i1, j, k;
+   int h, i, i1, j, k, nres;
    double agedeb;
    double ***p3mat;
- Line 9200  int hPijx(double *p, int bage, int fage)
+ Line 9379  int hPijx(double *p, int bage, int fage)
    /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
    /*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
    /*    k=k+1;  */
-   for (k=1; k <= (int) pow(2,cptcoveff); k++){
+   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-     fprintf(ficrespijb,"\n#****** ");
+     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-     for(j=1;j<=cptcoveff;j++)
+       if(TKresult[nres]!= k)
-       fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+         continue;
-     fprintf(ficrespijb,"******\n");
+       fprintf(ficrespijb,"\n#****** ");
-     if(invalidvarcomb[k]){
+       for(j=1;j<=cptcoveff;j++)
-       fprintf(ficrespijb,"\n#Combination (%d) ignored because no cases \n",k);
+         fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-       continue;
+       for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
-     }
+         fprintf(ficrespijb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+       }
-     /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */
+       fprintf(ficrespijb,"******\n");
-     for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */
+       if(invalidvarcomb[k]){
-       /* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /\* Typically 20 years = 20*12/6=40 *\/ */
+         fprintf(ficrespijb,"\n#Combination (%d) ignored because no cases \n",k);
-       nhstepm=(int) rint((agedeb-ageminl)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
+         continue;
-       nhstepm = nhstepm/hstepm; /* Typically 40/4=10, because estepm=24 stepm=6 => hstepm=24/6=4 */
+       }
-       /*          nhstepm=nhstepm*YEARM; aff par mois*/
+       /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */
+       for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */
-       p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+         /* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /\* Typically 20 years = 20*12/6=40 *\/ */
-       /* oldm=oldms;savm=savms; */
+         nhstepm=(int) rint((agedeb-ageminl)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
-       /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);   */
+         nhstepm = nhstepm/hstepm; /* Typically 40/4=10, because estepm=24 stepm=6 => hstepm=24/6=4 */
-       hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k);
-       /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
+         /*        nhstepm=nhstepm*YEARM; aff par mois*/
-       fprintf(ficrespijb,"# Cov Agex agex-h hpijx with i,j=");
-       for(i=1; i<=nlstate;i++)
+         p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
-         for(j=1; j<=nlstate+ndeath;j++)
+         /* oldm=oldms;savm=savms; */
-           fprintf(ficrespijb," %1d-%1d",i,j);
+         /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);   */
-       fprintf(ficrespijb,"\n");
+         hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k);
-       for (h=0; h<=nhstepm; h++){
+         /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
-         /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
+         fprintf(ficrespijb,"# Cov Agex agex-h hpijx with i,j=");
-         fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb - h*hstepm/YEARM*stepm );
-         /* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */
          for(i=1; i<=nlstate;i++)
            for(j=1; j<=nlstate+ndeath;j++)
-             fprintf(ficrespijb," %.5f", p3mat[i][j][h]);
+             fprintf(ficrespijb," %1d-%1d",i,j);
          fprintf(ficrespijb,"\n");
-       }
+         for (h=0; h<=nhstepm; h++){
-       free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+           /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
-       fprintf(ficrespijb,"\n");
+           fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb - h*hstepm/YEARM*stepm );
-     }
+           /* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */
-     /*}*/
+           for(i=1; i<=nlstate;i++)
-   }
+             for(j=1; j<=nlstate+ndeath;j++)
+               fprintf(ficrespijb," %.5f", p3mat[i][j][h]);
+           fprintf(ficrespijb,"\n");
+         }
+         free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+         fprintf(ficrespijb,"\n");
+       } /* end age deb */
+     } /* end combination */
+   } /* end nres */
    return 0;
   } /*  hBijx */
- Line 10646  Please run with mle=-1 to get a correct
+ Line 10831  Please run with mle=-1 to get a correct
        ungetc(c,ficpar);
        fgets(line, MAXLINE, ficpar);
        fputs(line,stdout);
+       fputs(line,ficres);
        fputs(line,ficparo);
      }
      ungetc(c,ficpar);
- Line 10662  Please run with mle=-1 to get a correct
+ Line 10848  Please run with mle=-1 to get a correct
        fgets(line, MAXLINE, ficpar);
        fputs(line,stdout);
        fputs(line,ficparo);
+       fputs(line,ficres);
      }
      ungetc(c,ficpar);
- Line 10680  Please run with mle=-1 to get a correct
+ Line 10867  Please run with mle=-1 to get a correct
          fputs(line,stdout);
          fputs(line,ficparo);
          fputs(line,ficlog);
+         fputs(line,ficres);
          continue;
        }else
          break;
      }
+     if (!feof(ficpar))
      while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){
-       if (num_filled == 0)
+       if (num_filled == 0){
          resultline[0]='\0';
-       else if (num_filled != 1){
+       break;
+       } else if (num_filled != 1){
          printf("ERROR %d: result line should be at minimum 'result=' %s\n",num_filled, line);
        }
        nresult++; /* Sum of resultlines */
- Line 10698  Please run with mle=-1 to get a correct
+ Line 10888  Please run with mle=-1 to get a correct
          goto end;
        }
        decoderesult(resultline, nresult); /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */
+       fprintf(ficparo,"result: %s\n",resultline);
+       fprintf(ficres,"result: %s\n",resultline);
+       fprintf(ficlog,"result: %s\n",resultline);
        while(fgets(line, MAXLINE, ficpar)) {
          /* If line starts with a # it is a comment */
          if (line[0] == '#') {
            numlinepar++;
            fputs(line,stdout);
            fputs(line,ficparo);
+           fputs(line,ficres);
            fputs(line,ficlog);
            continue;
          }else
- Line 10713  Please run with mle=-1 to get a correct
+ Line 10907  Please run with mle=-1 to get a correct
          break;
        else{ /* Processess output results for this combination of covariate values */
        }
-     }
+     } /* end while */
- Line 10848  Please run with mle=-1 to get a correct
+ Line 11042  Please run with mle=-1 to get a correct
      }
      printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout);
      fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog);
+     pstamp(ficreseij);
      i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
      if (cptcovn < 1){i1=1;}
- Line 10960  Please run with mle=-1 to get a correct
+ Line 11156  Please run with mle=-1 to get a correct
        fprintf(ficrescveij,"******\n");
        fprintf(ficresvij,"\n#****** ");
+       /* pstamp(ficresvij); */
        for(j=1;j<=cptcoveff;j++)
          fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
        for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
- Line 11137  Please run with mle=-1 to get a correct
+ Line 11334  Please run with mle=-1 to get a correct
    free_ivector(ncodemaxwundef,1,NCOVMAX);
    free_ivector(Dummy,-1,NCOVMAX);
    free_ivector(Fixed,-1,NCOVMAX);
+   free_ivector(DummyV,1,NCOVMAX);
+   free_ivector(FixedV,1,NCOVMAX);
    free_ivector(Typevar,-1,NCOVMAX);
    free_ivector(Tvar,1,NCOVMAX);
    free_ivector(TvarsQ,1,NCOVMAX);

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

Removed from v.1.236
changed lines
	Added in v.1.240