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

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

version 1.236, 2016/08/25 10:50:18 version 1.247, 2016/09/02 11:11:21
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.247  2016/09/02 11:11:21  brouard
     *** empty log message ***
   
     Revision 1.246  2016/09/02 08:49:22  brouard
     *** empty log message ***
   
     Revision 1.245  2016/09/02 07:25:01  brouard
     *** empty log message ***
   
     Revision 1.244  2016/09/02 07:17:34  brouard
     *** empty log message ***
   
     Revision 1.243  2016/09/02 06:45:35  brouard
     *** empty log message ***
   
     Revision 1.242  2016/08/30 15:01:20  brouard
     Summary: Fixing a lots
   
     Revision 1.241  2016/08/29 17:17:25  brouard
     Summary: gnuplot problem in Back projection to fix
   
     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    Revision 1.236  2016/08/25 10:50:18  brouard
   *** empty log message ***    *** empty log message ***
   
Line 1026  double dval; Line 1061  double dval;
 #define FTOL 1.0e-10  #define FTOL 1.0e-10
   
 #define NRANSI   #define NRANSI 
 #define ITMAX 200   #define ITMAX 200
   #define ITPOWMAX 20 /* This is now multiplied by the number of parameters */ 
   
 #define TOL 2.0e-4   #define TOL 2.0e-4 
   
Line 1114  int *TvarsQind; Line 1150  int *TvarsQind;
 #define MAXRESULTLINES 10  #define MAXRESULTLINES 10
 int nresult=0;  int nresult=0;
 int TKresult[MAXRESULTLINES];  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) */  int Tvresult[MAXRESULTLINES][NCOVMAX]; /* For dummy variable , variable # (output) */
 double Tqresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , value (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 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 *\/ */  /* 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 1178  double *Tvalsel; /**< Selected modality
 int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product */  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 *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 *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 *Tage;
 int anyvaryingduminmodel=0; /**< Any varying dummy in Model=1 yes, 0 no, to avoid a loop on waves in freq */   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*/  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 1189  int *Ndum; /** Freq of modality (tricode
 /* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */  /* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */
 int **Tvard;  int **Tvard;
 int *Tprod;/**< Gives the k position of the k1 product */  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 */  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;  int cptcovprod, *Tvaraff, *invalidvarcomb;
 double *lsurv, *lpop, *tpop;  double *lsurv, *lpop, *tpop;
   
Line 2052  void powell(double p[], double **xi, int Line 2091  void powell(double p[], double **xi, int
  void linmin(double p[], double xi[], int n, double *fret,    void linmin(double p[], double xi[], int n, double *fret, 
               double (*func)(double []));                 double (*func)(double [])); 
 #else   #else 
  void linmin(double p[], double xi[], int n, double *fret,    void linmin(double p[], double xi[], int n, double *fret,
                                                  double (*func)(double []),int *flat);                double (*func)(double []),int *flat); 
 #endif  #endif
   int i,ibig,j;    int i,ibig,j,jk,k; 
   double del,t,*pt,*ptt,*xit;    double del,t,*pt,*ptt,*xit;
   double directest;    double directest;
   double fp,fptt;    double fp,fptt;
Line 2087  void powell(double p[], double **xi, int Line 2126  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(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); */  /*     fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */
     for (i=1;i<=n;i++) {      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," %.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  ");
           fprintf(ficlog,"  + age*age  ");
       }
       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("\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(ficlog,"\n");
     fprintf(ficrespow,"\n");fflush(ficrespow);      for(i=1,jk=1; i <=nlstate; i++){
     if(*iter <=3){        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 && *iter >1){
       tml = *localtime(&rcurr_time);        tml = *localtime(&rcurr_time);
       strcpy(strcurr,asctime(&tml));        strcpy(strcurr,asctime(&tml));
       rforecast_time=rcurr_time;         rforecast_time=rcurr_time; 
       itmp = strlen(strcurr);        itmp = strlen(strcurr);
       if(strcurr[itmp-1]=='\n')  /* Windows outputs with a new line */        if(strcurr[itmp-1]=='\n')  /* Windows outputs with a new line */
                                 strcurr[itmp-1]='\0';          strcurr[itmp-1]='\0';
       printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);        printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
       fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);        fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
       for(niterf=10;niterf<=30;niterf+=10){        for(niterf=10;niterf<=30;niterf+=10){
                                 rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);          rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);
                                 forecast_time = *localtime(&rforecast_time);          forecast_time = *localtime(&rforecast_time);
                                 strcpy(strfor,asctime(&forecast_time));          strcpy(strfor,asctime(&forecast_time));
                                 itmp = strlen(strfor);          itmp = strlen(strfor);
                                 if(strfor[itmp-1]=='\n')          if(strfor[itmp-1]=='\n')
                                         strfor[itmp-1]='\0';            strfor[itmp-1]='\0';
                                 printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);          printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
                                 fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);          fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
       }        }
     }      }
     for (i=1;i<=n;i++) { /* For each direction i */      for (i=1;i<=n;i++) { /* For each direction i */
Line 2165  void powell(double p[], double **xi, int Line 2237  void powell(double p[], double **xi, int
                                 /* printf("\n"); */                                  /* printf("\n"); */
                                 /* fprintf(ficlog,"\n"); */                                  /* fprintf(ficlog,"\n"); */
                         }                          }
     if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */      /* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */
       if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */
       /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */        /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */
       /* By adding age*age in a model, the new -2LL should be lower and the difference follows a */        /* By adding age*age in a model, the new -2LL should be lower and the difference follows a */
       /* a chisquare statistics with 1 degree. To be significant at the 95% level, it should have */        /* a chisquare statistics with 1 degree. To be significant at the 95% level, it should have */
Line 2211  void powell(double p[], double **xi, int Line 2284  void powell(double p[], double **xi, int
       free_vector(pt,1,n);         free_vector(pt,1,n); 
       return;         return; 
     } /* enough precision */       } /* 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) */      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];         ptt[j]=2.0*p[j]-pt[j]; 
       xit[j]=p[j]-pt[j];         xit[j]=p[j]-pt[j]; 
Line 2410  void powell(double p[], double **xi, int Line 2483  void powell(double p[], double **xi, int
       cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];         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]); */        /* 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]]]){        if(Dummy[Tvar[Tage[k]]]){
         cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];          cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
       } else{        } else{
Line 2418  void powell(double p[], double **xi, int Line 2491  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]); */        /* 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]); */        /* 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 cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
     /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/      /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
Line 2479  Earliest age to start was %d-%d=%d, ncvl Line 2564  Earliest age to start was %d-%d=%d, ncvl
   
  /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ageminpar, double agemaxpar, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */   /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ageminpar, double agemaxpar, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */
  /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */   /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */
  double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ftolpl, int *ncvyear, int ij)    double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ftolpl, int *ncvyear, int ij, int nres)
 {  {
   /* Computes the prevalence limit in each live state at age x and covariate ij by left multiplying the unit    /* Computes the prevalence limit in each live state at age x and covariate ij by left multiplying the unit
      matrix by transitions matrix until convergence is reached with precision ftolpl */       matrix by transitions matrix until convergence is reached with precision ftolpl */
Line 2502  Earliest age to start was %d-%d=%d, ncvl Line 2587  Earliest age to start was %d-%d=%d, ncvl
   /* If we start from prlim again, prlim tends to a constant matrix */    /* If we start from prlim again, prlim tends to a constant matrix */
   
   int i, ii,j,k;    int i, ii,j,k;
     int first=0;
   double *min, *max, *meandiff, maxmax,sumnew=0.;    double *min, *max, *meandiff, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */    /* double **matprod2(); */ /* test */
   double **out, cov[NCOVMAX+1], **bmij();    double **out, cov[NCOVMAX+1], **bmij();
Line 2538  Earliest age to start was %d-%d=%d, ncvl Line 2624  Earliest age to start was %d-%d=%d, ncvl
     cov[2]=agefin;      cov[2]=agefin;
     if(nagesqr==1)      if(nagesqr==1)
       cov[3]= agefin*agefin;;        cov[3]= agefin*agefin;;
     for (k=1; k<=cptcovn;k++) {      for (k=1; k<=nsd;k++) { /* For single dummy covariates only */
       /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */                          /* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */
       cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];        cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];
       /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */        /* printf("bprevalim Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */
       }
       /* for (k=1; k<=cptcovn;k++) { */
       /*   /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */
       /*   cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */
       /*   /\* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); *\/ */
       /* } */
       for (k=1; k<=nsq;k++) { /* For single varying covariates only */
                           /* Here comes the value of quantitative after renumbering k with single quantitative covariates */
         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++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2]; */
       /* for (k=1; k<=cptcovprod;k++) /\* Useless *\/ */
       /*   /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; *\/ */
       /*   cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,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{
           cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; 
         }
         /* 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 product without age */
         /* printf("prevalim Prod ij=%d k=%d  Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */
         if(Dummy[Tvard[k][1]==0]){
           if(Dummy[Tvard[k][2]==0]){
             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]];
           }
         }
     }      }
     for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2];  
     for (k=1; k<=cptcovprod;k++) /* Useless */  
       /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */  
       cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];  
           
     /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/      /*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]);*/      /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
Line 2593  Earliest age to start was %d-%d=%d, ncvl Line 2713  Earliest age to start was %d-%d=%d, ncvl
     }      }
   } /* age loop */    } /* age loop */
     /* After some age loop it doesn't converge */      /* After some age loop it doesn't converge */
   printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\    if(first){
       first=1;
       printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. Others in log file only...\n\
   Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
     }
     fprintf(ficlog,"Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
 Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);  Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
   /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */    /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
   free_vector(min,1,nlstate);    free_vector(min,1,nlstate);
Line 2726  double **pmij(double **ps, double *cov, Line 2851  double **pmij(double **ps, double *cov,
         /* }else */          /* }else */
         doldm[ii][j]=(ii==j ? 1./sumnew : 0.0);          doldm[ii][j]=(ii==j ? 1./sumnew : 0.0);
       }else{        }else{
         printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin);          ;
           /* printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin); */
       }        }
     } /*End ii */      } /*End ii */
   } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */    } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */
Line 3109  double func( double *x) Line 3235  double func( double *x)
          Then computes with function pmij which return a matrix p[i][j] giving the elementary probability           Then computes with function pmij which return a matrix p[i][j] giving the elementary probability
          to be observed in j being in i according to the model.           to be observed in j being in i according to the model.
       */        */
       ioffset=2+nagesqr+cptcovage;        ioffset=2+nagesqr ;
    /* Fixed */     /* Fixed */
       for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */        for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */
         cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/          cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/
Line 3128  double func( double *x) Line 3254  double func( double *x)
       */        */
       for(mi=1; mi<= wav[i]-1; mi++){        for(mi=1; mi<= wav[i]-1; mi++){
         for(k=1; k <= ncovv ; k++){ /* Varying  covariates (single and product but no age )*/          for(k=1; k <= ncovv ; k++){ /* Varying  covariates (single and product but no age )*/
           cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i];            /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */
             cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
         }          }
         for (ii=1;ii<=nlstate+ndeath;ii++)          for (ii=1;ii<=nlstate+ndeath;ii++)
           for (j=1;j<=nlstate+ndeath;j++){            for (j=1;j<=nlstate+ndeath;j++){
Line 3142  double func( double *x) Line 3269  double func( double *x)
           if(nagesqr==1)            if(nagesqr==1)
             cov[3]= agexact*agexact;  /* Should be changed here */              cov[3]= agexact*agexact;  /* Should be changed here */
           for (kk=1; kk<=cptcovage;kk++) {            for (kk=1; kk<=cptcovage;kk++) {
             if(!FixedV[Tvar[Tage[kk]]])
             cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */              cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
             else
               cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact;
           }            }
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,            out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                        1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));                         1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
Line 3424  double funcone( double *x) Line 3554  double funcone( double *x)
   for(k=1; k<=nlstate; k++) ll[k]=0.;    for(k=1; k<=nlstate; k++) ll[k]=0.;
   ioffset=0;    ioffset=0;
   for (i=1,ipmx=0, sw=0.; i<=imx; i++){    for (i=1,ipmx=0, sw=0.; i<=imx; i++){
     ioffset=2+nagesqr+cptcovage;      /* ioffset=2+nagesqr+cptcovage; */
       ioffset=2+nagesqr;
     /* Fixed */      /* Fixed */
     /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */      /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */
     /* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */      /* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */
Line 3451  double funcone( double *x) Line 3582  double funcone( double *x)
     for(mi=1; mi<= wav[i]-1; mi++){  /* Varying with waves */      for(mi=1; mi<= wav[i]-1; mi++){  /* Varying with waves */
     /* Wave varying (but not age varying) */      /* Wave varying (but not age varying) */
       for(k=1; k <= ncovv ; k++){ /* Varying  covariates (single and product but no age )*/        for(k=1; k <= ncovv ; k++){ /* Varying  covariates (single and product but no age )*/
                                 cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i];          /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */
                         }          cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
         }
       /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */        /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */
                                 /* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */        /* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */
                                 /* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */        /* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */
                                 /* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */        /* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */
                                 /* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */        /* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */
                                 /* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */        /* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */
       /* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates *\/ */        /* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates *\/ */
                         /*      iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */        /*        iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */
                         /*      /\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */        /*        /\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */
                         /*      cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; */        /*        cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; */
       /* } */        /* } */
       for (ii=1;ii<=nlstate+ndeath;ii++)        for (ii=1;ii<=nlstate+ndeath;ii++)
                                 for (j=1;j<=nlstate+ndeath;j++){          for (j=1;j<=nlstate+ndeath;j++){
                                         oldm[ii][j]=(ii==j ? 1.0 : 0.0);            oldm[ii][j]=(ii==j ? 1.0 : 0.0);
                                         savm[ii][j]=(ii==j ? 1.0 : 0.0);            savm[ii][j]=(ii==j ? 1.0 : 0.0);
                                 }          }
               
       agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */        agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */
       ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */        ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */
       for(d=0; d<dh[mi][i]; d++){  /* Delay between two effective waves */        for(d=0; d<dh[mi][i]; d++){  /* Delay between two effective waves */
                                 /*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]        /* for(d=0; d<=0; d++){  /\* Delay between two effective waves Only one matrix to speed up*\/ */
                                         and mw[mi+1][i]. dh depends on stepm.*/          /*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
                                 newm=savm;            and mw[mi+1][i]. dh depends on stepm.*/
                                 agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;          newm=savm;
                                 cov[2]=agexact;          agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;  /* Here d is needed */
                                 if(nagesqr==1)          cov[2]=agexact;
                                         cov[3]= agexact*agexact;          if(nagesqr==1)
                                 for (kk=1; kk<=cptcovage;kk++) {            cov[3]= agexact*agexact;
                                         cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;          for (kk=1; kk<=cptcovage;kk++) {
                                 }            if(!FixedV[Tvar[Tage[kk]]])
                                 /* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */              cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
                                 /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */            else
                                 out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,              cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact;
                                                                                  1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));          }
                                 /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, */          /* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */
                                 /*           1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); */          /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
                                 savm=oldm;          out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                                 oldm=newm;                       1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
           /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, */
           /*           1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); */
           savm=oldm;
           oldm=newm;
       } /* end mult */        } /* end mult */
               
       s1=s[mw[mi][i]][i];        s1=s[mw[mi][i]][i];
Line 3504  double funcone( double *x) Line 3640  double funcone( double *x)
        * is higher than the multiple of stepm and negative otherwise.         * is higher than the multiple of stepm and negative otherwise.
        */         */
       if( s2 > nlstate && (mle <5) ){  /* Jackson */        if( s2 > nlstate && (mle <5) ){  /* Jackson */
                                 lli=log(out[s1][s2] - savm[s1][s2]);          lli=log(out[s1][s2] - savm[s1][s2]);
       } else if  ( s2==-1 ) { /* alive */        } else if  ( s2==-1 ) { /* alive */
                                 for (j=1,survp=0. ; j<=nlstate; j++)           for (j=1,survp=0. ; j<=nlstate; j++) 
                                         survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];            survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
                                 lli= log(survp);          lli= log(survp);
       }else if (mle==1){        }else if (mle==1){
                                 lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */          lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
       } else if(mle==2){        } else if(mle==2){
                                 lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */          lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */
       } else if(mle==3){  /* exponential inter-extrapolation */        } else if(mle==3){  /* exponential inter-extrapolation */
                                 lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */          lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
       } else if (mle==4){  /* mle=4 no inter-extrapolation */        } else if (mle==4){  /* mle=4 no inter-extrapolation */
                                 lli=log(out[s1][s2]); /* Original formula */          lli=log(out[s1][s2]); /* Original formula */
       } else{  /* mle=0 back to 1 */        } else{  /* mle=0 back to 1 */
                                 lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */          lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
                                 /*lli=log(out[s1][s2]); */ /* Original formula */          /*lli=log(out[s1][s2]); */ /* Original formula */
       } /* End of if */        } /* End of if */
       ipmx +=1;        ipmx +=1;
       sw += weight[i];        sw += weight[i];
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;        ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
       /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */        /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
       if(globpr){        if(globpr){
                                 fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\          fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
  %11.6f %11.6f %11.6f ", \   %11.6f %11.6f %11.6f ", \
                                                                 num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,                  num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
                                                                 2*weight[i]*lli,out[s1][s2],savm[s1][s2]);                  2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
                                 for(k=1,llt=0.,l=0.; k<=nlstate; k++){          for(k=1,llt=0.,l=0.; k<=nlstate; k++){
                                         llt +=ll[k]*gipmx/gsw;            llt +=ll[k]*gipmx/gsw;
                                         fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);            fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
                                 }          }
                                 fprintf(ficresilk," %10.6f\n", -llt);          fprintf(ficresilk," %10.6f\n", -llt);
       }        }
         } /* end of wave */          } /* end of wave */
 } /* end of individual */  } /* end of individual */
Line 3884  double hessij( double x[], double **hess Line 4020  double hessij( double x[], double **hess
       kmax=kmax+10;        kmax=kmax+10;
     }      }
     if(kmax >=10 || firstime ==1){      if(kmax >=10 || firstime ==1){
       printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol);        printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
       fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol);        fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
       printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);        printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
       fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);        fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
     }      }
Line 4073  void  freqsummary(char fileres[], int ia Line 4209  void  freqsummary(char fileres[], int ia
     fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);      fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
     exit(0);      exit(0);
   }    }
     
   strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));    strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));
   if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {    if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {
     printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));      printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
Line 4083  void  freqsummary(char fileres[], int ia Line 4219  void  freqsummary(char fileres[], int ia
   }    }
   else{    else{
     fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \      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",\  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);              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"));    strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
   if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {    if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
     printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));      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));      fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
     fflush(ficlog);      fflush(ficlog);
     exit(70);       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> \      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",\  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);              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);    freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
   j1=0;    j1=0;
       
   /* j=ncoveff;  /\* Only fixed dummy covariates *\/ */    /* j=ncoveff;  /\* Only fixed dummy covariates *\/ */
   j=cptcoveff;  /* Only dummy covariates of the model */    j=cptcoveff;  /* Only dummy covariates of the model */
   if (cptcovn<1) {j=1;ncodemax[1]=1;}    if (cptcovn<1) {j=1;ncodemax[1]=1;}
     
   first=1;    first=1;
     
   /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels:    /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels:
      reference=low_education V1=0,V2=0       reference=low_education V1=0,V2=0
      med_educ                V1=1 V2=0,        med_educ                V1=1 V2=0, 
      high_educ               V1=0 V2=1       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        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 */    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.;      posproptt=0.;
     /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);      /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
       scanf("%d", i);*/        scanf("%d", i);*/
     for (i=-5; i<=nlstate+ndeath; i++)        for (i=-5; i<=nlstate+ndeath; i++)  
       for (jk=-5; jk<=nlstate+ndeath; jk++)          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 (i=1; i<=nlstate; i++)  {
       for(m=iagemin; m <= iagemax+3; m++)        for(m=iagemin; m <= iagemax+3; m++)
                                 prop[i][m]=0;          prop[i][m]=0;
       posprop[i]=0;        posprop[i]=0;
       pospropt[i]=0;        pospropt[i]=0;
     }      }
Line 4141  Title=%s <br>Datafile=%s Firstpass=%d La Line 4276  Title=%s <br>Datafile=%s Firstpass=%d La
     /*  meanqt[m][z1]=0.; */      /*  meanqt[m][z1]=0.; */
     /*   } */      /*   } */
     /* } */      /* } */
                       
     dateintsum=0;      dateintsum=0;
     k2cpt=0;      k2cpt=0;
     /* For that combination of covariate j1, we count and print the frequencies in one pass */      /* 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 4355  Title=%s <br>Datafile=%s Firstpass=%d La
     } /* end iind = 1 to imx */      } /* end iind = 1 to imx */
     /* prop[s][age] is feeded for any initial and valid live state as well as      /* 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 */         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);*/      /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
     pstamp(ficresp);      pstamp(ficresp);
     /* if  (ncoveff>0) { */      if  (cptcoveff>0){
     if  (cptcoveff>0) {  
       fprintf(ficresp, "\n#********** Variable ");         fprintf(ficresp, "\n#********** Variable "); 
       fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");         fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); 
       fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");         fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); 
         fprintf(ficlog, "\n#********** Variable "); 
       for (z1=1; z1<=cptcoveff; z1++){        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(ficresp, "**********\n#");
       fprintf(ficresphtm, "**********</h3>\n");        fprintf(ficresphtm, "**********</h3>\n");
       fprintf(ficresphtmfr, "**********</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(ficlog, "**********\n");
     }      }
     fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");      fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
     for(i=1; i<=nlstate;i++) {      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(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);
     }      }
     fprintf(ficresp, "\n");      fprintf(ficresp, "\n");
     fprintf(ficresphtm, "\n");      fprintf(ficresphtm, "\n");
                       
     /* Header of frequency table by age */      /* Header of frequency table by age */
     fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");      fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");
     fprintf(ficresphtmfr,"<th>Age</th> ");      fprintf(ficresphtmfr,"<th>Age</th> ");
Line 4259  Title=%s <br>Datafile=%s Firstpass=%d La Line 4400  Title=%s <br>Datafile=%s Firstpass=%d La
       }        }
     }      }
     fprintf(ficresphtmfr, "\n");      fprintf(ficresphtmfr, "\n");
                       
     /* For each age */      /* For each age */
     for(iage=iagemin; iage <= iagemax+3; iage++){      for(iage=iagemin; iage <= iagemax+3; iage++){
       fprintf(ficresphtm,"<tr>");        fprintf(ficresphtm,"<tr>");
       if(iage==iagemax+1){        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){        }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){        }else if(iage==iagemax+3){
                                 fprintf(ficlog,"Total");          fprintf(ficlog,"Total");
                                 fprintf(ficresphtmfr,"<tr><th>Total</th> ");          fprintf(ficresphtmfr,"<tr><th>Total</th> ");
       }else{        }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(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(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++){         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 */        } /* 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++){        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++){        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 */        } /* end loop jk */
       /* pospropt=0.; */        /* pospropt=0.; */
       for(jk=-1; jk <=nlstate+ndeath; jk++){        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 */        } /* end loop jk */
       posproptt=0.;         posproptt=0.; 
       for(jk=1; jk <=nlstate; jk++){        for(jk=1; jk <=nlstate; jk++){
                                 posproptt += pospropt[jk];          posproptt += pospropt[jk];
       }        }
       fprintf(ficresphtmfr,"</tr>\n ");        fprintf(ficresphtmfr,"</tr>\n ");
       if(iage <= iagemax){        if(iage <= iagemax){
                                 fprintf(ficresp,"\n");          fprintf(ficresp,"\n");
                                 fprintf(ficresphtm,"</tr>\n");          fprintf(ficresphtm,"</tr>\n");
       }        }
       if(first==1)        if(first==1)
                                 printf("Others in log...\n");          printf("Others in log...\n");
       fprintf(ficlog,"\n");        fprintf(ficlog,"\n");
     } /* end loop age iage */      } /* end loop age iage */
     fprintf(ficresphtm,"<tr><th>Tot</th>");      fprintf(ficresphtm,"<tr><th>Tot</th>");
     for(jk=1; jk <=nlstate ; jk++){      for(jk=1; jk <=nlstate ; jk++){
       if(posproptt < 1.e-5){        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{        }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");      fprintf(ficresphtm,"</tr>\n");
Line 4385  Title=%s <br>Datafile=%s Firstpass=%d La Line 4526  Title=%s <br>Datafile=%s Firstpass=%d La
     fprintf(ficresphtmfr,"</table>\n");      fprintf(ficresphtmfr,"</table>\n");
   } /* end selected combination of covariate j1 */    } /* end selected combination of covariate j1 */
   dateintmean=dateintsum/k2cpt;     dateintmean=dateintsum/k2cpt; 
             
   fclose(ficresp);    fclose(ficresp);
   fclose(ficresphtm);    fclose(ficresphtm);
   fclose(ficresphtmfr);    fclose(ficresphtmfr);
Line 4720  void  concatwav(int wav[], int **dh, int Line 4861  void  concatwav(int wav[], int **dh, int
   
 /*********** Tricode ****************************/  /*********** Tricode ****************************/
  void tricode(int *cptcov, int *Tvar, int **nbcode, int imx, int *Ndum)   void tricode(int *cptcov, int *Tvar, int **nbcode, int imx, int *Ndum)
 {   {
   /**< Uses cptcovn+2*cptcovprod as the number of covariates */     /**< Uses cptcovn+2*cptcovprod as the number of covariates */
   /*      Tvar[i]=atoi(stre);  find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1      /*     Tvar[i]=atoi(stre);  find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 
    * Boring subroutine which should only output nbcode[Tvar[j]][k]      * Boring subroutine which should only output nbcode[Tvar[j]][k]
    * Tvar[5] in V2+V1+V3*age+V2*V4 is 4 (V4) even it is a time varying or quantitative variable      * Tvar[5] in V2+V1+V3*age+V2*V4 is 4 (V4) even it is a time varying or quantitative variable
    * nbcode[Tvar[5]][1]= nbcode[4][1]=0, nbcode[4][2]=1 (usually);      * nbcode[Tvar[5]][1]= nbcode[4][1]=0, nbcode[4][2]=1 (usually);
   */      */
   
   int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;     int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
   int modmaxcovj=0; /* Modality max of covariates j */     int modmaxcovj=0; /* Modality max of covariates j */
   int cptcode=0; /* Modality max of covariates j */     int cptcode=0; /* Modality max of covariates j */
   int modmincovj=0; /* Modality min of covariates j */     int modmincovj=0; /* Modality min of covariates j */
   
   
   /* cptcoveff=0;  */     /* cptcoveff=0;  */
         /* *cptcov=0; */     /* *cptcov=0; */
     
   for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */     for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */
   
   /* Loop on covariates without age and products and no quantitative variable */     /* Loop on covariates without age and products and no quantitative variable */
   /* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */     /* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */
   for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */     for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */
     for (j=-1; (j < maxncov); j++) Ndum[j]=0;       for (j=-1; (j < maxncov); j++) Ndum[j]=0;
     if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */        if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ 
       switch(Fixed[k]) {         switch(Fixed[k]) {
       case 0: /* Testing on fixed dummy covariate, simple or product of fixed */         case 0: /* Testing on fixed dummy covariate, simple or product of fixed */
                                 for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the  modality of this covariate Vj*/           for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the  modality of this covariate Vj*/
                                         ij=(int)(covar[Tvar[k]][i]);             ij=(int)(covar[Tvar[k]][i]);
                                         /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i             /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
                                          * If product of Vn*Vm, still boolean *:              * If product of Vn*Vm, still boolean *:
                                          * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables              * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
                                          * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0   */              * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0   */
                                         /* Finds for covariate j, n=Tvar[j] of Vn . ij is the             /* Finds for covariate j, n=Tvar[j] of Vn . ij is the
                                                  modality of the nth covariate of individual i. */                modality of the nth covariate of individual i. */
                                         if (ij > modmaxcovj)             if (ij > modmaxcovj)
                                                 modmaxcovj=ij;                modmaxcovj=ij; 
                                         else if (ij < modmincovj)              else if (ij < modmincovj) 
                                                 modmincovj=ij;                modmincovj=ij; 
                                         if ((ij < -1) && (ij > NCOVMAX)){             if ((ij < -1) && (ij > NCOVMAX)){
                                                 printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );               printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
                                                 exit(1);               exit(1);
                                         }else             }else
                                                 Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/               Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/
                                         /*  If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */             /*  If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */
                                         /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/             /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
                                         /* getting the maximum value of the modality of the covariate             /* getting the maximum value of the modality of the covariate
                                                  (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and                (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and
                                                  female ies 1, then modmaxcovj=1.                female ies 1, then modmaxcovj=1.
                                         */             */
                                 } /* end for loop on individuals i */           } /* end for loop on individuals i */
                                 printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);           printf(" Minimal and maximal values of %d th (fixed) covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);
                                 fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);           fprintf(ficlog," Minimal and maximal values of %d th (fixed) covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);
                                 cptcode=modmaxcovj;           cptcode=modmaxcovj;
                                 /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */           /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */
                                 /*for (i=0; i<=cptcode; i++) {*/           /*for (i=0; i<=cptcode; i++) {*/
                                 for (j=modmincovj;  j<=modmaxcovj; j++) { /* j=-1 ? 0 and 1*//* For each value j of the modality of model-cov k */           for (j=modmincovj;  j<=modmaxcovj; j++) { /* j=-1 ? 0 and 1*//* For each value j of the modality of model-cov k */
                                         printf("Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);             printf("Frequencies of (fixed) covariate %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);
                                         fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);             fprintf(ficlog, "Frequencies of (fixed) covariate %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);
                                         if( Ndum[j] != 0 ){ /* Counts if nobody answered modality j ie empty modality, we skip it and reorder */             if( Ndum[j] != 0 ){ /* Counts if nobody answered modality j ie empty modality, we skip it and reorder */
                                                 if( j != -1){               if( j != -1){
                                                         ncodemax[k]++;  /* ncodemax[k]= Number of modalities of the k th                 ncodemax[k]++;  /* ncodemax[k]= Number of modalities of the k th
                                                                                                                                  covariate for which somebody answered excluding                                     covariate for which somebody answered excluding 
                                                                                                                                  undefined. Usually 2: 0 and 1. */                                    undefined. Usually 2: 0 and 1. */
                                                 }               }
                                                 ncodemaxwundef[k]++; /* ncodemax[j]= Number of modalities of the k th               ncodemaxwundef[k]++; /* ncodemax[j]= Number of modalities of the k th
                                                                                                                                                 covariate for which somebody answered including                                        covariate for which somebody answered including 
                                                                                                                                                 undefined. Usually 3: -1, 0 and 1. */                                       undefined. Usually 3: -1, 0 and 1. */
                                         }       /* In fact  ncodemax[k]=2 (dichotom. variables only) but it could be more for             }    /* In fact  ncodemax[k]=2 (dichotom. variables only) but it could be more for
                                                  * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */                   * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
                                 } /* Ndum[-1] number of undefined modalities */           } /* Ndum[-1] number of undefined modalities */
                                                   
                                 /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */           /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
                                 /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. */           /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. */
                                 /* If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125; */           /* If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125; */
                                 /* modmincovj=3; modmaxcovj = 7; */           /* modmincovj=3; modmaxcovj = 7; */
                                 /* There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3; */           /* There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3; */
                                 /* which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10; */           /* which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10; */
                           /*             defining two dummy variables: variables V1_1 and V1_2.*/           /*              defining two dummy variables: variables V1_1 and V1_2.*/
               /* nbcode[Tvar[j]][ij]=k; */           /* nbcode[Tvar[j]][ij]=k; */
               /* nbcode[Tvar[j]][1]=0; */           /* nbcode[Tvar[j]][1]=0; */
               /* nbcode[Tvar[j]][2]=1; */           /* nbcode[Tvar[j]][2]=1; */
               /* nbcode[Tvar[j]][3]=2; */           /* nbcode[Tvar[j]][3]=2; */
               /* To be continued (not working yet). */           /* To be continued (not working yet). */
               ij=0; /* ij is similar to i but can jump over null modalities */           ij=0; /* ij is similar to i but can jump over null modalities */
                                 for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/           for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/
           if (Ndum[i] == 0) { /* If nobody responded to this modality k */             if (Ndum[i] == 0) { /* If nobody responded to this modality k */
                   break;               break;
                 }             }
                                         ij++;             ij++;
                                         nbcode[Tvar[k]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/             nbcode[Tvar[k]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/
                                         cptcode = ij; /* New max modality for covar j */             cptcode = ij; /* New max modality for covar j */
                                 } /* end of loop on modality i=-1 to 1 or more */           } /* end of loop on modality i=-1 to 1 or more */
                                 break;           break;
       case 1: /* Testing on varying covariate, could be simple and         case 1: /* Testing on varying covariate, could be simple and
                * should look at waves or product of fixed *                  * should look at waves or product of fixed *
                * varying. No time to test -1, assuming 0 and 1 only */                  * varying. No time to test -1, assuming 0 and 1 only */
                                 ij=0;           ij=0;
                                 for(i=0; i<=1;i++){           for(i=0; i<=1;i++){
                                         nbcode[Tvar[k]][++ij]=i;             nbcode[Tvar[k]][++ij]=i;
                                 }           }
                                 break;           break;
       default:         default:
                                 break;           break;
       } /* end switch */         } /* end switch */
     } /* end dummy test */       } /* end dummy test */
           
     /*   for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */       /*   for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */
     /*  /\*recode from 0 *\/ */       /*         /\*recode from 0 *\/ */
     /*                               k is a modality. If we have model=V1+V1*sex  */       /*                                      k is a modality. If we have model=V1+V1*sex  */
     /*                               then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */       /*                                      then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */
     /*                            But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */       /*                                   But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */
     /*  } */       /*         } */
     /*  /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */       /*         /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */
     /*  if (ij > ncodemax[j]) { */       /*         if (ij > ncodemax[j]) { */
     /*    printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]);  */       /*           printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]);  */
     /*    fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */       /*           fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */
     /*    break; */       /*           break; */
     /*  } */       /*         } */
     /*   }  /\* end of loop on modality k *\/ */       /*   }  /\* end of loop on modality k *\/ */
   } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/       } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/  
       
   for (k=-1; k< maxncov; k++) Ndum[k]=0;      for (k=-1; k< maxncov; k++) Ndum[k]=0; 
   /* Look at fixed dummy (single or product) covariates to check empty modalities */     /* Look at fixed dummy (single or product) covariates to check empty modalities */
   for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */      for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */ 
     /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/        /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ 
     ij=Tvar[i]; /* Tvar 5,4,3,6,5,7,1,4 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V4*age */        ij=Tvar[i]; /* Tvar 5,4,3,6,5,7,1,4 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V4*age */ 
     Ndum[ij]++; /* Count the # of 1, 2 etc: {1,1,1,2,2,1,1} because V1 once, V2 once, two V4 and V5 in above */       Ndum[ij]++; /* Count the # of 1, 2 etc: {1,1,1,2,2,1,1} because V1 once, V2 once, two V4 and V5 in above */
     /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1,  {2, 1, 1, 1, 2, 1, 1, 0, 0} */       /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1,  {2, 1, 1, 1, 2, 1, 1, 0, 0} */
   } /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 1} */     } /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 1} */
       
   ij=0;     ij=0;
   /* for (i=0; i<=  maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) *\/ */     /* for (i=0; i<=  maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) *\/ */
   for (k=1; k<=  cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */     for (k=1; k<=  cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
     /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/       /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
     /* if((Ndum[i]!=0) && (i<=ncovcol)){  /\* Tvar[i] <= ncovmodel ? *\/ */       /* if((Ndum[i]!=0) && (i<=ncovcol)){  /\* Tvar[i] <= ncovmodel ? *\/ */
     if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){  /* Only Dummy and non empty in the model */       if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){  /* Only Dummy and non empty in the model */
       /* If product not in single variable we don't print results */         /* If product not in single variable we don't print results */
       /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/         /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
       ++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */         ++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */
       Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/         Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/
       Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */         Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */
       TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */         TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */
       if(Fixed[k]!=0)         if(Fixed[k]!=0)
         anyvaryingduminmodel=1;           anyvaryingduminmodel=1;
                         /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */         /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */
                         /*   Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet *\/ */         /*   Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet *\/ */
                         /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */         /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */
                         /*   Tvaraff[++ij]=i; /\*For printing (unclear) *\/ */         /*   Tvaraff[++ij]=i; /\*For printing (unclear) *\/ */
                         /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */         /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */
                         /*   Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet *\/ */         /*   Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet *\/ */
     }        } 
   } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */     } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */
   /* ij--; */     /* ij--; */
   /* cptcoveff=ij; /\*Number of total covariates*\/ */     /* cptcoveff=ij; /\*Number of total covariates*\/ */
   *cptcov=ij; /*Number of total real effective covariates: effective     *cptcov=ij; /*Number of total real effective covariates: effective
                                                          * because they can be excluded from the model and real                  * because they can be excluded from the model and real
                                                          * if in the model but excluded because missing values, but how to get k from ij?*/                  * if in the model but excluded because missing values, but how to get k from ij?*/
   for(j=ij+1; j<= cptcovt; j++){     for(j=ij+1; j<= cptcovt; j++){
     Tvaraff[j]=0;       Tvaraff[j]=0;
     Tmodelind[j]=0;       Tmodelind[j]=0;
   }     }
   for(j=ntveff+1; j<= cptcovt; j++){     for(j=ntveff+1; j<= cptcovt; j++){
     TmodelInvind[j]=0;       TmodelInvind[j]=0;
   }     }
   /* To be sorted */     /* To be sorted */
   ;     ;
 }   }
   
   
 /*********** Health Expectancies ****************/  /*********** Health Expectancies ****************/
Line 4901  void  concatwav(int wav[], int **dh, int Line 5042  void  concatwav(int wav[], int **dh, int
   double ***p3mat;    double ***p3mat;
   double eip;    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,"# (a) Life expectancies by health status at initial age and (b) health expectancies by health status at initial age\n");
   fprintf(ficreseij,"# Age");    fprintf(ficreseij,"# Age");
   for(i=1; i<=nlstate;i++){    for(i=1; i<=nlstate;i++){
Line 5266  void  concatwav(int wav[], int **dh, int Line 5407  void  concatwav(int wav[], int **dh, int
    fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);     fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
    pstamp(ficresprobmorprev);     pstamp(ficresprobmorprev);
    fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);     fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
      fprintf(ficresprobmorprev,"# Selected quantitative variables and dummies");
      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);     fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
    for(j=nlstate+1; j<=(nlstate+ndeath);j++){     for(j=nlstate+1; j<=(nlstate+ndeath);j++){
      fprintf(ficresprobmorprev," p.%-d SE",j);       fprintf(ficresprobmorprev," p.%-d SE",j);
Line 5335  void  concatwav(int wav[], int **dh, int Line 5484  void  concatwav(int wav[], int **dh, int
          xp[i] = x[i] + (i==theta ?delti[theta]:0);           xp[i] = x[i] + (i==theta ?delti[theta]:0);
        }         }
                                                   
        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nresult);         prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres);
                                                   
        if (popbased==1) {         if (popbased==1) {
          if(mobilav ==0){           if(mobilav ==0){
Line 5367  void  concatwav(int wav[], int **dh, int Line 5516  void  concatwav(int wav[], int **dh, int
        for(i=1; i<=npar; i++) /* Computes gradient x - delta */         for(i=1; i<=npar; i++) /* Computes gradient x - delta */
          xp[i] = x[i] - (i==theta ?delti[theta]:0);           xp[i] = x[i] - (i==theta ?delti[theta]:0);
                                                   
        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nresult);         prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nres);
                                                   
        if (popbased==1) {         if (popbased==1) {
          if(mobilav ==0){           if(mobilav ==0){
Line 5444  void  concatwav(int wav[], int **dh, int Line 5593  void  concatwav(int wav[], int **dh, int
      /* end ppptj */       /* end ppptj */
      /*  x centered again */       /*  x centered again */
                                   
      prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij, nresult);       prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij, nres);
                                   
      if (popbased==1) {       if (popbased==1) {
        if(mobilav ==0){         if(mobilav ==0){
Line 5539  void  concatwav(int wav[], int **dh, int Line 5688  void  concatwav(int wav[], int **dh, int
       
   pstamp(ficresvpl);    pstamp(ficresvpl);
   fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n");    fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n");
   fprintf(ficresvpl,"# Age");    fprintf(ficresvpl,"# Age ");
     if(nresult >=1)
       fprintf(ficresvpl," Result# ");
   for(i=1; i<=nlstate;i++)    for(i=1; i<=nlstate;i++)
       fprintf(ficresvpl," %1d-%1d",i,i);        fprintf(ficresvpl," %1d-%1d",i,i);
   fprintf(ficresvpl,"\n");    fprintf(ficresvpl,"\n");
Line 5625  void  concatwav(int wav[], int **dh, int Line 5776  void  concatwav(int wav[], int **dh, int
       varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */        varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */
   
     fprintf(ficresvpl,"%.0f ",age );      fprintf(ficresvpl,"%.0f ",age );
       if(nresult >=1)
         fprintf(ficresvpl,"%d ",nres );
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
       fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age]));        fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age]));
     fprintf(ficresvpl,"\n");      fprintf(ficresvpl,"\n");
Line 5986  void printinghtml(char fileresu[], char Line 6139  void printinghtml(char fileresu[], char
                   int popforecast, int prevfcast, int backcast, int estepm , \                    int popforecast, int prevfcast, int backcast, int estepm , \
                   double jprev1, double mprev1,double anprev1, double dateprev1, \                    double jprev1, double mprev1,double anprev1, double dateprev1, \
                   double jprev2, double mprev2,double anprev2, double dateprev2){                    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 \     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 \     <li><a href='#secondorder'>Result files (second order (variance)</a>\n \
 </ul>");  </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,"<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",     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"));             jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm"));
Line 6025  void printinghtml(char fileresu[], char Line 6180  void printinghtml(char fileresu[], char
    if (cptcovn < 1) {m=1;ncodemax[1]=1;}     if (cptcovn < 1) {m=1;ncodemax[1]=1;}
   
    jj1=0;     jj1=0;
    for(k1=1; k1<=m;k1++){  
      for(nres=1; nres <= nresult; nres++) /* For each resultline */
      for(k1=1; k1<=m;k1++){ /* For each combination of covariate */
        if(TKresult[nres]!= k1)
          continue;
   
      /* for(i1=1; i1<=ncodemax[k1];i1++){ */       /* for(i1=1; i1<=ncodemax[k1];i1++){ */
      jj1++;       jj1++;
      if (cptcovn > 0) {       if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");         fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++){          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 */         /* if(nqfveff+nqtveff 0) */ /* Test to be done */
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
        if(invalidvarcomb[k1]){         if(invalidvarcomb[k1]){
Line 6044  void printinghtml(char fileresu[], char Line 6210  void printinghtml(char fileresu[], char
        }         }
      }       }
      /* aij, bij */       /* aij, bij */
      fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \       fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1-%d.svg\">%s_%d-1-%d.svg</a><br> \
 <img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);  <img src=\"%s_%d-1-%d.svg\">",model,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
      /* Pij */       /* Pij */
      fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \       fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2-%d.svg\">%s_%d-2-%d.svg</a><br> \
 <img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);       <img src=\"%s_%d-2-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);     
      /* Quasi-incidences */       /* Quasi-incidences */
      fprintf(fichtm,"<br>\n- I<sub>ij</sub> or Conditional probabilities to be observed in state j being in state i %d (stepm) months\       fprintf(fichtm,"<br>\n- I<sub>ij</sub> or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
  before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too, \   before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too, \
  incidence (rates) are the limit when h tends to zero of the ratio of the probability  <sub>h</sub>P<sub>ij</sub> \   incidence (rates) are the limit when h tends to zero of the ratio of the probability  <sub>h</sub>P<sub>ij</sub> \
 divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \  divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3-%d.svg\">%s_%d-3-%d.svg</a><br> \
 <img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);   <img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); 
      /* Survival functions (period) in state j */       /* Survival functions (period) in state j */
      for(cpt=1; cpt<=nlstate;cpt++){       for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \         fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
 <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1);  <img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
      }       }
      /* State specific survival functions (period) */       /* State specific survival functions (period) */
      for(cpt=1; cpt<=nlstate;cpt++){       for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\         fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\
  Or probability to survive in various states (1 to %d) being in state %d at different ages.     \   Or probability to survive in various states (1 to %d) being in state %d at different ages.     \
  <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> <img src=\"%s_%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1);   <a href=\"%s_%d-%d-%d.svg\">%s_%d%d-%d.svg</a><br> <img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
      }       }
      /* Period (stable) prevalence in each health state */       /* Period (stable) prevalence in each health state */
      for(cpt=1; cpt<=nlstate;cpt++){       for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \         fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
 <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);  <img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
      }       }
      if(backcast==1){       if(backcast==1){
        /* Period (stable) back prevalence in each health state */         /* Period (stable) back prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){         for(cpt=1; cpt<=nlstate;cpt++){
          fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \           fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
 <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1);  <img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
        }         }
      }       }
      if(prevfcast==1){       if(prevfcast==1){
        /* Projection of prevalence up to period (stable) prevalence in each health state */         /* Projection of prevalence up to period (stable) prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){         for(cpt=1; cpt<=nlstate;cpt++){
          fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \           fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
 <img src=\"%s_%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1);  <img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
        }         }
      }       }
                     
      for(cpt=1; cpt<=nlstate;cpt++) {       for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \         fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a> <br> \
 <img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1);  <img src=\"%s_%d-%d-%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres);
      }       }
      /* } /\* end i1 *\/ */       /* } /\* end i1 *\/ */
    }/* End k1 */     }/* End k1 */
Line 6146  See page 'Matrix of variance-covariance Line 6312  See page 'Matrix of variance-covariance
    if (cptcovn < 1) {m=1;ncodemax[1]=1;}     if (cptcovn < 1) {m=1;ncodemax[1]=1;}
   
    jj1=0;     jj1=0;
   
      for(nres=1; nres <= nresult; nres++){ /* For each resultline */
    for(k1=1; k1<=m;k1++){     for(k1=1; k1<=m;k1++){
        if(TKresult[nres]!= k1)
          continue;
      /* for(i1=1; i1<=ncodemax[k1];i1++){ */       /* for(i1=1; i1<=ncodemax[k1];i1++){ */
      jj1++;       jj1++;
      if (cptcovn > 0) {       if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");         fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++)  /**< cptcoveff number of variables */         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\">");         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
   
        if(invalidvarcomb[k1]){         if(invalidvarcomb[k1]){
Line 6162  See page 'Matrix of variance-covariance Line 6337  See page 'Matrix of variance-covariance
      }       }
      for(cpt=1; cpt<=nlstate;cpt++) {       for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \         fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \
 prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d.svg\"> %s_%d-%d.svg</a>\n <br>\  prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d-%d.svg\"> %s_%d-%d-%d.svg</a>\n <br>\
 <img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1);    <img src=\"%s_%d-%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres);  
      }       }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and \       fprintf(fichtm,"\n<br>- Total life expectancy by age and \
 health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \  health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
 true period expectancies (those weighted with period prevalences are also\  true period expectancies (those weighted with period prevalences are also\
  drawn in addition to the population based expectancies computed using\   drawn in addition to the population based expectancies computed using\
  observed and cahotic prevalences:  <a href=\"%s_%d.svg\">%s_%d.svg</a>\n<br>\   observed and cahotic prevalences:  <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>\n<br>\
 <img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1);  <img src=\"%s_%d-%d.svg\">",subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres);
      /* } /\* end i1 *\/ */       /* } /\* end i1 *\/ */
    }/* End k1 */     }/* End k1 */
     }/* End nres */
    fprintf(fichtm,"</ul>");     fprintf(fichtm,"</ul>");
    fflush(fichtm);     fflush(fichtm);
 }  }
Line 6182  void printinggnuplot(char fileresu[], ch Line 6358  void printinggnuplot(char fileresu[], ch
   
   char dirfileres[132],optfileres[132];    char dirfileres[132],optfileres[132];
   char gplotcondition[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 lv=0, vlv=0, kl=0;
   int ng=0;    int ng=0;
   int vpopbased;    int vpopbased;
Line 6231  void printinggnuplot(char fileresu[], ch Line 6407  void printinggnuplot(char fileresu[], ch
   strcpy(dirfileres,optionfilefiname);    strcpy(dirfileres,optionfilefiname);
   strcpy(optfileres,"vpl");    strcpy(optfileres,"vpl");
   /* 1eme*/    /* 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 */        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-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
       fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);          fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);
       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==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres);
 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==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres);
         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==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); 
         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, nres in 2 to be fixed */
       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 lt 3",    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], \  
                       4+(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++;
             }                /* 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' w l lt 3",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
                           2+cptcoveff*2+(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 */      } /* k1 */
   } /* cpt */    } /* cpt */
   
       
   /*2 eme*/    /*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 */        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 */          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,1,4) = 1 because h=1  k= (1) 1  1  1 */
Line 6390  plot [%.f:%.f] \"%s\" every :::%d::%d u Line 6504  plot [%.f:%.f] \"%s\" every :::%d::%d u
         vlv= nbcode[Tvaraff[k]][lv];          vlv= nbcode[Tvaraff[k]][lv];
         fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);          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 */        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," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
       }         }
       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;
       }        }
                                                   
       /*       k=2+nlstate*(2*cpt-2); */        fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1,nres);
       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);
 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);          else
       /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);            fprintf(ficgp,"\nreplot ");
         for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");          for (i=1; i<= nlstate+1 ; i ++) {
         fprintf(ficgp,"\" t \"e%d1\" w l",cpt);            k=2*i;
         fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);            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 (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");            for (j=1; j<= nlstate+1 ; j ++) {
         fprintf(ficgp,"\" t \"e%d1\" w l",cpt);              if (j==i) fprintf(ficgp," %%lf (%%lf)");
                                               else fprintf(ficgp," %%*lf (%%*lf)");
       */            }   
       for (i=1; i< nlstate ; i ++) {            if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",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);            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:%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:($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)");
       fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);              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);
   /* 4eme */            for (j=1; j<= nlstate+1 ; j ++) {
   /* Survival functions (period) from state i in state j by initial state i */              if (j==i) fprintf(ficgp," %%lf (%%lf)");
   for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */              else fprintf(ficgp," %%*lf (%%*lf)");
   for(nres=1; nres <= nresult; nres++) /* For each resultline */            }   
     if(TKresult[nres]!= k)            if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
       continue;            else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
           } /* state */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */        } /* vpopbased */
       fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);        fprintf(ficgp,"\nset out;set out \"%s_%d-%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */
         lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */      } /* end nres */
         /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */    } /* k1 end 2 eme*/
         /* 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];    /*3eme*/
         fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);    for (k1=1; k1<= m ; k1 ++){
       }      for(nres=1; nres <= nresult; nres++){ /* For each resultline */
       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */        if(TKresult[nres]!= k1)
         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;          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;
           }
                                                   
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);          /*       k=2+nlstate*(2*cpt-2); */
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\          k=2+(nlstate+1)*(cpt-1);
 set ter svg size 640, 480\n                                                                                                                                                                                     \          fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres);
 unset log y\n                                                                                                                                                                                                                                           \          fprintf(ficgp,"set ter svg size 640, 480\n\
 plot [%.f:%.f]  ", ageminpar, agemaxpar);  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);
       k=3;          /*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 ; i ++){            for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
         if(i==1){            fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
           fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));            fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
         }else{            for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
           fprintf(ficgp,", '' ");            fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
                                   
           */
           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(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);
         }
       } /* end nres */
     } /* end kl 3eme */
     
     /* 4eme */
     /* Survival functions (period) from state i in state j by initial state i */
     for (k1=1; k1<=m; k1++){    /* For each covariate and each value */
       for(nres=1; nres <= nresult; nres++){ /* For each resultline */
         if(TKresult[nres]!= 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);
           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);
         }          }
         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-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
           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 */  /* 5eme */
   /* Survival functions (period) from state i in state j by final state j */    /* 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;          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-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
       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 */  /* 6eme */
   /* CV preval stable (period) for each covariate */    /* 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 */      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);        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 6715  plot [%.f:%.f]  ", ageminpar, agemaxpar)
         vlv= nbcode[Tvaraff[k]][lv];          vlv= nbcode[Tvaraff[k]][lv];
         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 */
           fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
         } 
       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,"P_"),cpt,k1);        fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\        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 */        k=3; /* Offset */
       for (i=1; i<= nlstate ; i ++){        for (i=1; i<= nlstate ; i ++){
         if(i==1)          if(i==1)
Line 6569  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6747  plot [%.f:%.f]  ", ageminpar, agemaxpar)
 /* 7eme */  /* 7eme */
   if(backcast == 1){    if(backcast == 1){
     /* CV back preval stable (period) for each covariate */      /* 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 */        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);          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 */          for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
Line 6580  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6761  plot [%.f:%.f]  ", ageminpar, agemaxpar)
           vlv= nbcode[Tvaraff[k]][lv];            vlv= nbcode[Tvaraff[k]][lv];
           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 */
             fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
           }       
         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,"PB_"),cpt,k1);          fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\          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 */          k=3; /* Offset */
         for (i=1; i<= nlstate ; i ++){          for (i=1; i<= nlstate ; i ++){
           if(i==1)            if(i==1)
Line 6616  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6798  plot [%.f:%.f]  ", ageminpar, agemaxpar)
   if(prevfcast==1){    if(prevfcast==1){
     /* Projection from cross-sectional to stable (period) for each covariate */      /* 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 */        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);          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  */          for (k=1; k<=cptcoveff; k++){    /* For each correspondig covariate value  */
Line 6627  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6812  plot [%.f:%.f]  ", ageminpar, agemaxpar)
           vlv= nbcode[Tvaraff[k]][lv];            vlv= nbcode[Tvaraff[k]][lv];
           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 */
             fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
           }       
         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); 
Line 6634  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6822  plot [%.f:%.f]  ", ageminpar, agemaxpar)
         }          }
                   
         fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");          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,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\          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 */          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*/            /*#  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 */               /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */   
Line 6705  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6891  plot [%.f:%.f]  ", ageminpar, agemaxpar)
   } /* End if prevfcast */    } /* 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++){    for(i=1,jk=1; i <=nlstate; i++){
     fprintf(ficgp,"# initial state %d\n",i);      fprintf(ficgp,"# initial state %d\n",i);
     for(k=1; k <=(nlstate+ndeath); k++){      for(k=1; k <=(nlstate+ndeath); k++){
Line 6723  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6909  plot [%.f:%.f]  ", ageminpar, agemaxpar)
   fprintf(ficgp,"##############\n#\n");    fprintf(ficgp,"##############\n#\n");
       
   /*goto avoid;*/    /*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)=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(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");    fprintf(ficgp,"# logi(p13/p11)=a13+b13*age+c13age*age+d13*V1+e13*V1*age\n");
Line 6738  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6925  plot [%.f:%.f]  ", ageminpar, agemaxpar)
   fprintf(ficgp,"#       +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n");    fprintf(ficgp,"#       +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n");
   fprintf(ficgp,"#\n");    fprintf(ficgp,"#\n");
   for(ng=1; ng<=3;ng++){ /* Number of graphics: first is logit, 2nd is probabilities, third is incidences per year*/    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 */
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng);      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-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng,nres);
       fprintf(ficgp,"\nset ter svg size 640, 480 ");        fprintf(ficgp,"\nset ter svg size 640, 480 ");
       if (ng==1){        if (ng==1){
         fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */          fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */
Line 6782  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6978  plot [%.f:%.f]  ", ageminpar, agemaxpar)
               break;                break;
             }              }
             ij=1;/* To be checked else nbcode[0][0] wrong */              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); */                /* 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++;                    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{            }else{
             i=i-ncovmodel;              i=i-ncovmodel;
             if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */              if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
Line 6811  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 7039  plot [%.f:%.f]  ", ageminpar, agemaxpar)
                                 
               ij=1;                ij=1;
               for(j=3; j <=ncovmodel-nagesqr; j++){                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,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])]); */                      /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
                     ij++;                      ij++;
Line 7076  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 7304  plot [%.f:%.f]  ", ageminpar, agemaxpar)
       fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);        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 */      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");      fprintf(ficresf," yearproj age");
     for(j=1; j<=nlstate+ndeath;j++){       for(j=1; j<=nlstate+ndeath;j++){ 
Line 7686  int readdata(char datafile[], int firsto Line 7913  int readdata(char datafile[], int firsto
   /*-------- data file ----------*/    /*-------- data file ----------*/
   FILE *fic;    FILE *fic;
   char dummy[]="                         ";    char dummy[]="                         ";
   int i=0, j=0, n=0, iv=0;    int i=0, j=0, n=0, iv=0, v;
   int lstra;    int lstra;
   int linei, month, year,iout;    int linei, month, year,iout;
   char line[MAXLINE], linetmp[MAXLINE];    char line[MAXLINE], linetmp[MAXLINE];
   char stra[MAXLINE], strb[MAXLINE];    char stra[MAXLINE], strb[MAXLINE];
   char *stratrunc;    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)    {    if((fic=fopen(datafile,"r"))==NULL)    {
     printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);      printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);
Line 7723  int readdata(char datafile[], int firsto Line 7972  int readdata(char datafile[], int firsto
     /* Loops on waves */      /* Loops on waves */
     for (j=maxwav;j>=1;j--){      for (j=maxwav;j>=1;j--){
       for (iv=nqtv;iv>=1;iv--){  /* Loop  on time varying quantitative variables */        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 */        }/* end loop ntqv */
               
       for (iv=ntv;iv>=1;iv--){  /* Loop  on time varying dummies */        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 \   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 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 \   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);   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 \   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 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 \   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);   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 */        }/* end loop ntv */
               
       /* Statuses  at wave */        /* Statuses  at wave */
       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 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;        s[j][i]=lval;
Line 8020  int decoderesult ( char resultline[], in Line 8269  int decoderesult ( char resultline[], in
     if(Typevar[k1]==0){ /* Single covariate in model */      if(Typevar[k1]==0){ /* Single covariate in model */
       match=0;        match=0;
       for(k2=1; k2 <=j;k2++){/* result line V4=1 V5=24.1 V3=1  V2=8 V1=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 */            modelresult[k2]=k1;/* modelresult[2]=1 modelresult[1]=2  modelresult[3]=3  modelresult[6]=4 modelresult[9]=5 */
           match=1;            match=1;
           break;            break;
Line 8036  int decoderesult ( char resultline[], in Line 8285  int decoderesult ( char resultline[], in
     match=0;      match=0;
     for(k1=1; k1<= cptcovt ;k1++){ /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */      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(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 */            resultmodel[k1]=k2;  /* resultmodel[2]=1 resultmodel[1]=2  resultmodel[3]=3  resultmodel[6]=4 resultmodel[9]=5 */
           ++match;            ++match;
         }          }
Line 8063  int decoderesult ( char resultline[], in Line 8312  int decoderesult ( char resultline[], in
   /*    6 1 0 1 */ /* V4=1, V3=0, V1=1 */    /*    6 1 0 1 */ /* V4=1, V3=0, V1=1 */
   /*    7 0 1 1 */    /*    7 0 1 1 */
   /*    8 1 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 */    for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* model line */
     if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Single dummy */      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 */        k2=(int)Tvarsel[k3]; /*  Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */
       k+=Tvalsel[k3]*pow(2,k4);  /*  Tvalsel[1]=1  */        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);        printf("Decoderesult Dummy k=%d, V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k2, k3, (int)Tvalsel[k3], k4);
       k4++;;        k4++;;
     }  else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Single quantitative */      }  else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Single quantitative */
       k3q= resultmodel[k1]; /* resultmodel[2] = 1=k3 */        k3q= resultmodel[k1]; /* resultmodel[2] = 1=k3 */
       k2q=(int)Tvarsel[k3q]; /*  Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */        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]);        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++;;        k4q++;;
     }      }
Line 8102  int decodemodel( char model[], int lasto Line 8357  int decodemodel( char model[], int lasto
         * - Tvard[k]  p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 .          * - 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;    int  j1, k1, k2, k3, k4;
   char modelsav[80];    char modelsav[80];
   char stra[80], strb[80], strc[80], strd[80],stre[80];    char stra[80], strb[80], strc[80], strd[80],stre[80];
Line 8308  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8563  Dummy[k] 0=dummy (0 1), 1 quantitative (
 Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product \n\  Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product \n\
 Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\  Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
 Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);  Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
     for(k=1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;}
   for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */    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 */      if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
       Fixed[k]= 0;        Fixed[k]= 0;
Line 8333  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8588  Dummy[k] 0=dummy (0 1), 1 quantitative (
       TvarFind[ncovf]=k;        TvarFind[ncovf]=k;
       TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */        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 */        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;        Fixed[k]= 0;
       Dummy[k]= 1;        Dummy[k]= 1;
       nqfveff++;        nqfveff++;
Line 8347  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8602  Dummy[k] 0=dummy (0 1), 1 quantitative (
       TvarFind[ncovf]=k;        TvarFind[ncovf]=k;
       TvarFQ[nqfveff]=Tvar[k]-ncovcol; /* TvarFQ[1]=V2-1=1st in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */        TvarFQ[nqfveff]=Tvar[k]-ncovcol; /* TvarFQ[1]=V2-1=1st in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */
       TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */        TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */
     }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying variables */      }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying dummy variables */
       Fixed[k]= 1;        Fixed[k]= 1;
       Dummy[k]= 0;        Dummy[k]= 0;
       ntveff++; /* Only simple time varying dummy variable */        ntveff++; /* Only simple time varying dummy variable */
Line 8358  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8613  Dummy[k] 0=dummy (0 1), 1 quantitative (
       TvarsDind[nsd]=k;        TvarsDind[nsd]=k;
       ncovv++; /* Only simple time varying variables */        ncovv++; /* Only simple time varying variables */
       TvarV[ncovv]=Tvar[k];        TvarV[ncovv]=Tvar[k];
       TvarVind[ncovv]=k;        TvarVind[ncovv]=k; /* TvarVind[2]=2  TvarVind[3]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */
       TvarVD[ntveff]=Tvar[k]; /* TvarVD[1]=V4  TvarVD[2]=V3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */        TvarVD[ntveff]=Tvar[k]; /* TvarVD[1]=V4  TvarVD[2]=V3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */
       TvarVDind[ntveff]=k; /* TvarVDind[1]=2 TvarVDind[2]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */        TvarVDind[ntveff]=k; /* TvarVDind[1]=2 TvarVDind[2]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying dummy variable */
       printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv);        printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv);
Line 8374  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8629  Dummy[k] 0=dummy (0 1), 1 quantitative (
       TvarsQ[nsq]=Tvar[k];        TvarsQ[nsq]=Tvar[k];
       TvarsQind[nsq]=k;        TvarsQind[nsq]=k;
       TvarV[ncovv]=Tvar[k];        TvarV[ncovv]=Tvar[k];
       TvarVind[ncovv]=k;        TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */
       TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */        TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */
       TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */        TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */
       TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */        TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */
Line 8386  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8641  Dummy[k] 0=dummy (0 1), 1 quantitative (
       TvarA[ncova]=Tvar[k];        TvarA[ncova]=Tvar[k];
       TvarAind[ncova]=k;        TvarAind[ncova]=k;
       if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */        if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */
         Fixed[k]= 2;          Fixed[k]= 2;
         Dummy[k]= 2;          Dummy[k]= 2;
         modell[k].maintype= ATYPE;          modell[k].maintype= ATYPE;
         modell[k].subtype= APFD;          modell[k].subtype= APFD;
         /* ncoveff++; */          /* ncoveff++; */
       }else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/        }else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/
         Fixed[k]= 2;          Fixed[k]= 2;
         Dummy[k]= 3;          Dummy[k]= 3;
         modell[k].maintype= ATYPE;          modell[k].maintype= ATYPE;
         modell[k].subtype= APFQ;                /*      Product age * fixed quantitative */          modell[k].subtype= APFQ;                /*      Product age * fixed quantitative */
         /* nqfveff++;  /\* Only simple fixed quantitative variable *\/ */          /* nqfveff++;  /\* Only simple fixed quantitative variable *\/ */
       }else if( Tvar[k] <=ncovcol+nqv+ntv ){        }else if( Tvar[k] <=ncovcol+nqv+ntv ){
         Fixed[k]= 3;          Fixed[k]= 3;
         Dummy[k]= 2;          Dummy[k]= 2;
         modell[k].maintype= ATYPE;          modell[k].maintype= ATYPE;
         modell[k].subtype= APVD;                /*      Product age * varying dummy */          modell[k].subtype= APVD;                /*      Product age * varying dummy */
         /* ntveff++; /\* Only simple time varying dummy variable *\/ */          /* ntveff++; /\* Only simple time varying dummy variable *\/ */
       }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){        }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){
         Fixed[k]= 3;          Fixed[k]= 3;
         Dummy[k]= 3;          Dummy[k]= 3;
         modell[k].maintype= ATYPE;          modell[k].maintype= ATYPE;
         modell[k].subtype= APVQ;                /*      Product age * varying quantitative */          modell[k].subtype= APVQ;                /*      Product age * varying quantitative */
         /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */          /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */
       }        }
     }else if (Typevar[k] == 2) {  /* product without age */      }else if (Typevar[k] == 2) {  /* product without age */
       k1=Tposprod[k];        k1=Tposprod[k];
       if(Tvard[k1][1] <=ncovcol){        if(Tvard[k1][1] <=ncovcol){
         if(Tvard[k1][2] <=ncovcol){          if(Tvard[k1][2] <=ncovcol){
           Fixed[k]= 1;            Fixed[k]= 1;
           Dummy[k]= 0;            Dummy[k]= 0;
           modell[k].maintype= FTYPE;            modell[k].maintype= FTYPE;
           modell[k].subtype= FPDD;              /*      Product fixed dummy * fixed dummy */            modell[k].subtype= FPDD;              /*      Product fixed dummy * fixed dummy */
           ncovf++; /* Fixed variables without age */            ncovf++; /* Fixed variables without age */
           TvarF[ncovf]=Tvar[k];            TvarF[ncovf]=Tvar[k];
           TvarFind[ncovf]=k;            TvarFind[ncovf]=k;
         }else if(Tvard[k1][2] <=ncovcol+nqv){          }else if(Tvard[k1][2] <=ncovcol+nqv){
           Fixed[k]= 0;  /* or 2 ?*/            Fixed[k]= 0;  /* or 2 ?*/
           Dummy[k]= 1;            Dummy[k]= 1;
           modell[k].maintype= FTYPE;            modell[k].maintype= FTYPE;
           modell[k].subtype= FPDQ;              /*      Product fixed dummy * fixed quantitative */            modell[k].subtype= FPDQ;              /*      Product fixed dummy * fixed quantitative */
           ncovf++; /* Varying variables without age */            ncovf++; /* Varying variables without age */
           TvarF[ncovf]=Tvar[k];            TvarF[ncovf]=Tvar[k];
           TvarFind[ncovf]=k;            TvarFind[ncovf]=k;
         }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
           Fixed[k]= 1;            Fixed[k]= 1;
           Dummy[k]= 0;            Dummy[k]= 0;
           modell[k].maintype= VTYPE;            modell[k].maintype= VTYPE;
           modell[k].subtype= VPDD;              /*      Product fixed dummy * varying dummy */            modell[k].subtype= VPDD;              /*      Product fixed dummy * varying dummy */
           ncovv++; /* Varying variables without age */            ncovv++; /* Varying variables without age */
           TvarV[ncovv]=Tvar[k];            TvarV[ncovv]=Tvar[k];
           TvarVind[ncovv]=k;            TvarVind[ncovv]=k;
         }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
           Fixed[k]= 1;            Fixed[k]= 1;
           Dummy[k]= 1;            Dummy[k]= 1;
           modell[k].maintype= VTYPE;            modell[k].maintype= VTYPE;
           modell[k].subtype= VPDQ;              /*      Product fixed dummy * varying quantitative */            modell[k].subtype= VPDQ;              /*      Product fixed dummy * varying quantitative */
           ncovv++; /* Varying variables without age */            ncovv++; /* Varying variables without age */
           TvarV[ncovv]=Tvar[k];            TvarV[ncovv]=Tvar[k];
           TvarVind[ncovv]=k;            TvarVind[ncovv]=k;
         }           }
       }else if(Tvard[k1][1] <=ncovcol+nqv){        }else if(Tvard[k1][1] <=ncovcol+nqv){
         if(Tvard[k1][2] <=ncovcol){          if(Tvard[k1][2] <=ncovcol){
           Fixed[k]= 0;  /* or 2 ?*/            Fixed[k]= 0;  /* or 2 ?*/
           Dummy[k]= 1;            Dummy[k]= 1;
           modell[k].maintype= FTYPE;            modell[k].maintype= FTYPE;
           modell[k].subtype= FPDQ;              /*      Product fixed quantitative * fixed dummy */            modell[k].subtype= FPDQ;              /*      Product fixed quantitative * fixed dummy */
           ncovf++; /* Fixed variables without age */            ncovf++; /* Fixed variables without age */
           TvarF[ncovf]=Tvar[k];            TvarF[ncovf]=Tvar[k];
           TvarFind[ncovf]=k;            TvarFind[ncovf]=k;
         }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
           Fixed[k]= 1;            Fixed[k]= 1;
           Dummy[k]= 1;            Dummy[k]= 1;
           modell[k].maintype= VTYPE;            modell[k].maintype= VTYPE;
           modell[k].subtype= VPDQ;              /*      Product fixed quantitative * varying dummy */            modell[k].subtype= VPDQ;              /*      Product fixed quantitative * varying dummy */
           ncovv++; /* Varying variables without age */            ncovv++; /* Varying variables without age */
           TvarV[ncovv]=Tvar[k];            TvarV[ncovv]=Tvar[k];
           TvarVind[ncovv]=k;            TvarVind[ncovv]=k;
         }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
           Fixed[k]= 1;            Fixed[k]= 1;
           Dummy[k]= 1;            Dummy[k]= 1;
           modell[k].maintype= VTYPE;            modell[k].maintype= VTYPE;
           modell[k].subtype= VPQQ;              /*      Product fixed quantitative * varying quantitative */            modell[k].subtype= VPQQ;              /*      Product fixed quantitative * varying quantitative */
           ncovv++; /* Varying variables without age */            ncovv++; /* Varying variables without age */
           TvarV[ncovv]=Tvar[k];            TvarV[ncovv]=Tvar[k];
           TvarVind[ncovv]=k;            TvarVind[ncovv]=k;
           ncovv++; /* Varying variables without age */            ncovv++; /* Varying variables without age */
           TvarV[ncovv]=Tvar[k];            TvarV[ncovv]=Tvar[k];
           TvarVind[ncovv]=k;            TvarVind[ncovv]=k;
         }           }
       }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){        }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){
         if(Tvard[k1][2] <=ncovcol){          if(Tvard[k1][2] <=ncovcol){
           Fixed[k]= 1;            Fixed[k]= 1;
           Dummy[k]= 1;            Dummy[k]= 1;
           modell[k].maintype= VTYPE;            modell[k].maintype= VTYPE;
           modell[k].subtype= VPDD;              /*      Product time varying dummy * fixed dummy */            modell[k].subtype= VPDD;              /*      Product time varying dummy * fixed dummy */
           ncovv++; /* Varying variables without age */            ncovv++; /* Varying variables without age */
           TvarV[ncovv]=Tvar[k];            TvarV[ncovv]=Tvar[k];
           TvarVind[ncovv]=k;            TvarVind[ncovv]=k;
         }else if(Tvard[k1][2] <=ncovcol+nqv){          }else if(Tvard[k1][2] <=ncovcol+nqv){
           Fixed[k]= 1;            Fixed[k]= 1;
           Dummy[k]= 1;            Dummy[k]= 1;
           modell[k].maintype= VTYPE;            modell[k].maintype= VTYPE;
           modell[k].subtype= VPDQ;              /*      Product time varying dummy * fixed quantitative */            modell[k].subtype= VPDQ;              /*      Product time varying dummy * fixed quantitative */
           ncovv++; /* Varying variables without age */            ncovv++; /* Varying variables without age */
           TvarV[ncovv]=Tvar[k];            TvarV[ncovv]=Tvar[k];
           TvarVind[ncovv]=k;            TvarVind[ncovv]=k;
         }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
           Fixed[k]= 1;            Fixed[k]= 1;
           Dummy[k]= 0;            Dummy[k]= 0;
           modell[k].maintype= VTYPE;            modell[k].maintype= VTYPE;
           modell[k].subtype= VPDD;              /*      Product time varying dummy * time varying dummy */            modell[k].subtype= VPDD;              /*      Product time varying dummy * time varying dummy */
           ncovv++; /* Varying variables without age */            ncovv++; /* Varying variables without age */
           TvarV[ncovv]=Tvar[k];            TvarV[ncovv]=Tvar[k];
           TvarVind[ncovv]=k;            TvarVind[ncovv]=k;
         }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
           Fixed[k]= 1;            Fixed[k]= 1;
           Dummy[k]= 1;            Dummy[k]= 1;
           modell[k].maintype= VTYPE;            modell[k].maintype= VTYPE;
           modell[k].subtype= VPDQ;              /*      Product time varying dummy * time varying quantitative */            modell[k].subtype= VPDQ;              /*      Product time varying dummy * time varying quantitative */
           ncovv++; /* Varying variables without age */            ncovv++; /* Varying variables without age */
           TvarV[ncovv]=Tvar[k];            TvarV[ncovv]=Tvar[k];
           TvarVind[ncovv]=k;            TvarVind[ncovv]=k;
         }           }
       }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){        }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){
         if(Tvard[k1][2] <=ncovcol){          if(Tvard[k1][2] <=ncovcol){
           Fixed[k]= 1;            Fixed[k]= 1;
           Dummy[k]= 1;            Dummy[k]= 1;
           modell[k].maintype= VTYPE;            modell[k].maintype= VTYPE;
           modell[k].subtype= VPDQ;              /*      Product time varying quantitative * fixed dummy */            modell[k].subtype= VPDQ;              /*      Product time varying quantitative * fixed dummy */
           ncovv++; /* Varying variables without age */            ncovv++; /* Varying variables without age */
           TvarV[ncovv]=Tvar[k];            TvarV[ncovv]=Tvar[k];
           TvarVind[ncovv]=k;            TvarVind[ncovv]=k;
         }else if(Tvard[k1][2] <=ncovcol+nqv){          }else if(Tvard[k1][2] <=ncovcol+nqv){
           Fixed[k]= 1;            Fixed[k]= 1;
           Dummy[k]= 1;            Dummy[k]= 1;
           modell[k].maintype= VTYPE;            modell[k].maintype= VTYPE;
           modell[k].subtype= VPQQ;              /*      Product time varying quantitative * fixed quantitative */            modell[k].subtype= VPQQ;              /*      Product time varying quantitative * fixed quantitative */
           ncovv++; /* Varying variables without age */            ncovv++; /* Varying variables without age */
           TvarV[ncovv]=Tvar[k];            TvarV[ncovv]=Tvar[k];
           TvarVind[ncovv]=k;            TvarVind[ncovv]=k;
         }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
           Fixed[k]= 1;            Fixed[k]= 1;
           Dummy[k]= 1;            Dummy[k]= 1;
           modell[k].maintype= VTYPE;            modell[k].maintype= VTYPE;
           modell[k].subtype= VPDQ;              /*      Product time varying quantitative * time varying dummy */            modell[k].subtype= VPDQ;              /*      Product time varying quantitative * time varying dummy */
           ncovv++; /* Varying variables without age */            ncovv++; /* Varying variables without age */
           TvarV[ncovv]=Tvar[k];            TvarV[ncovv]=Tvar[k];
           TvarVind[ncovv]=k;            TvarVind[ncovv]=k;
         }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
           Fixed[k]= 1;            Fixed[k]= 1;
           Dummy[k]= 1;            Dummy[k]= 1;
           modell[k].maintype= VTYPE;            modell[k].maintype= VTYPE;
           modell[k].subtype= VPQQ;              /*      Product time varying quantitative * time varying quantitative */            modell[k].subtype= VPQQ;              /*      Product time varying quantitative * time varying quantitative */
           ncovv++; /* Varying variables without age */            ncovv++; /* Varying variables without age */
           TvarV[ncovv]=Tvar[k];            TvarV[ncovv]=Tvar[k];
           TvarVind[ncovv]=k;            TvarVind[ncovv]=k;
         }           }
       }else{        }else{
         printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);          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]);          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{      }else{
       printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);        printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
       fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);        fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
Line 8923  int prevalence_limit(double *p, double * Line 9178  int prevalence_limit(double *p, double *
   i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */    i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
   if (cptcovn < 1){i1=1;}    if (cptcovn < 1){i1=1;}
   
   for(nres=1; nres <= nresult; nres++) /* For each resultline */    for(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;    return 0;
 }  }
   
Line 9018  int back_prevalence_limit(double *p, dou Line 9275  int back_prevalence_limit(double *p, dou
   i1=pow(2,cptcoveff);    i1=pow(2,cptcoveff);
   if (cptcovn < 1){i1=1;}    if (cptcovn < 1){i1=1;}
       
   for(nres=1; nres <= nresult; nres++) /* For each resultline */    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
   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, nres);
       }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,nres);
       }          }
       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); */    /* hBijx(p, bage, fage); */
   /* fclose(ficrespijb); */    /* fclose(ficrespijb); */
       
Line 9172  int hPijx(double *p, int bage, int fage) Line 9429  int hPijx(double *p, int bage, int fage)
         int ageminl;          int ageminl;
   int hstepm;    int hstepm;
   int nhstepm;    int nhstepm;
   int h, i, i1, j, k;    int h, i, i1, j, k, nres;
                   
   double agedeb;    double agedeb;
   double ***p3mat;    double ***p3mat;
Line 9200  int hPijx(double *p, int bage, int fage) Line 9457  int hPijx(double *p, int bage, int fage)
   /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */    /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
   /*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */    /*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
   /*    k=k+1;  */    /*    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(i=1; i<=nlstate;i++)
           for(j=1; j<=nlstate+ndeath;j++)            for(j=1; j<=nlstate+ndeath;j++)
             fprintf(ficrespijb," %.5f", p3mat[i][j][h]);              fprintf(ficrespijb," %1d-%1d",i,j);
         fprintf(ficrespijb,"\n");          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;    return 0;
  } /*  hBijx */   } /*  hBijx */
   
Line 9624  int main(int argc, char *argv[]) Line 9887  int main(int argc, char *argv[])
   delti=delti3[1][1];    delti=delti3[1][1];
   /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/    /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/
   if(mle==-1){ /* Print a wizard for help writing covariance matrix */    if(mle==-1){ /* Print a wizard for help writing covariance matrix */
   /* We could also provide initial parameters values giving by simple logistic regression 
    * only one way, that is without matrix product. We will have nlstate maximizations */
         /* for(i=1;i<nlstate;i++){ */
         /*        /\*reducing xi for 1 to npar to 1 to ncovmodel; *\/ */
         /*    mlikeli(ficres,p, ncovmodel, ncovmodel, nlstate, ftol, funcnoprod); */
         /* } */
     prwizard(ncovmodel, nlstate, ndeath, model, ficparo);      prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
     printf(" You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);      printf(" You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
     fprintf(ficlog," You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);      fprintf(ficlog," You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
Line 10396  Please run with mle=-1 to get a correct Line 10665  Please run with mle=-1 to get a correct
     printf("\n");      printf("\n");
     if(mle>=1){ /* Could be 1 or 2, Real Maximization */      if(mle>=1){ /* Could be 1 or 2, Real Maximization */
       /* mlikeli uses func not funcone */        /* mlikeli uses func not funcone */
         /* for(i=1;i<nlstate;i++){ */
         /*        /\*reducing xi for 1 to npar to 1 to ncovmodel; *\/ */
         /*    mlikeli(ficres,p, ncovmodel, ncovmodel, nlstate, ftol, funcnoprod); */
         /* } */
       mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);        mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
     }      }
     if(mle==0) {/* No optimization, will print the likelihoods for the datafile */      if(mle==0) {/* No optimization, will print the likelihoods for the datafile */
Line 10646  Please run with mle=-1 to get a correct Line 10919  Please run with mle=-1 to get a correct
       ungetc(c,ficpar);        ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);        fgets(line, MAXLINE, ficpar);
       fputs(line,stdout);        fputs(line,stdout);
         fputs(line,ficres);
       fputs(line,ficparo);        fputs(line,ficparo);
     }      }
     ungetc(c,ficpar);      ungetc(c,ficpar);
Line 10662  Please run with mle=-1 to get a correct Line 10936  Please run with mle=-1 to get a correct
       fgets(line, MAXLINE, ficpar);        fgets(line, MAXLINE, ficpar);
       fputs(line,stdout);        fputs(line,stdout);
       fputs(line,ficparo);        fputs(line,ficparo);
         fputs(line,ficres);
     }      }
     ungetc(c,ficpar);      ungetc(c,ficpar);
           
Line 10680  Please run with mle=-1 to get a correct Line 10955  Please run with mle=-1 to get a correct
         fputs(line,stdout);          fputs(line,stdout);
         fputs(line,ficparo);          fputs(line,ficparo);
         fputs(line,ficlog);          fputs(line,ficlog);
           fputs(line,ficres);
         continue;          continue;
       }else        }else
         break;          break;
     }      }
       if (!feof(ficpar))
     while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){      while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){
       if (num_filled == 0)        if (num_filled == 0){
         resultline[0]='\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);          printf("ERROR %d: result line should be at minimum 'result=' %s\n",num_filled, line);
       }        }
       nresult++; /* Sum of resultlines */        nresult++; /* Sum of resultlines */
Line 10698  Please run with mle=-1 to get a correct Line 10976  Please run with mle=-1 to get a correct
         goto end;          goto end;
       }        }
       decoderesult(resultline, nresult); /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */        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)) {        while(fgets(line, MAXLINE, ficpar)) {
         /* If line starts with a # it is a comment */          /* If line starts with a # it is a comment */
         if (line[0] == '#') {          if (line[0] == '#') {
           numlinepar++;            numlinepar++;
           fputs(line,stdout);            fputs(line,stdout);
           fputs(line,ficparo);            fputs(line,ficparo);
             fputs(line,ficres);
           fputs(line,ficlog);            fputs(line,ficlog);
           continue;            continue;
         }else          }else
Line 10713  Please run with mle=-1 to get a correct Line 10995  Please run with mle=-1 to get a correct
         break;          break;
       else{ /* Processess output results for this combination of covariate values */        else{ /* Processess output results for this combination of covariate values */
       }                                    }                            
     }      } /* end while */
   
   
           
Line 10848  Please run with mle=-1 to get a correct Line 11130  Please run with mle=-1 to get a correct
     }      }
     printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout);      printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout);
     fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog);      fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog);
   
       pstamp(ficreseij);
                                   
     i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */      i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
     if (cptcovn < 1){i1=1;}      if (cptcovn < 1){i1=1;}
Line 10929  Please run with mle=-1 to get a correct Line 11213  Please run with mle=-1 to get a correct
     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("\n#****** Selected:");        printf("\n#****** Result for:");
       fprintf(ficrest,"\n#****** Selected:");        fprintf(ficrest,"\n#****** Result for:");
       fprintf(ficlog,"\n#****** Selected:");        fprintf(ficlog,"\n#****** Result for:");
       for(j=1;j<=cptcoveff;j++){         for(j=1;j<=cptcoveff;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(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);          fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
Line 10960  Please run with mle=-1 to get a correct Line 11244  Please run with mle=-1 to get a correct
       fprintf(ficrescveij,"******\n");        fprintf(ficrescveij,"******\n");
               
       fprintf(ficresvij,"\n#****** ");        fprintf(ficresvij,"\n#****** ");
         /* pstamp(ficresvij); */
       for(j=1;j<=cptcoveff;j++)         for(j=1;j<=cptcoveff;j++) 
         fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,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 */        for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
Line 11137  Please run with mle=-1 to get a correct Line 11422  Please run with mle=-1 to get a correct
   free_ivector(ncodemaxwundef,1,NCOVMAX);    free_ivector(ncodemaxwundef,1,NCOVMAX);
   free_ivector(Dummy,-1,NCOVMAX);    free_ivector(Dummy,-1,NCOVMAX);
   free_ivector(Fixed,-1,NCOVMAX);    free_ivector(Fixed,-1,NCOVMAX);
     free_ivector(DummyV,1,NCOVMAX);
     free_ivector(FixedV,1,NCOVMAX);
   free_ivector(Typevar,-1,NCOVMAX);    free_ivector(Typevar,-1,NCOVMAX);
   free_ivector(Tvar,1,NCOVMAX);    free_ivector(Tvar,1,NCOVMAX);
   free_ivector(TvarsQ,1,NCOVMAX);    free_ivector(TvarsQ,1,NCOVMAX);
Removed from v.1.236  
changed lines
  Added in v.1.247

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