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

Diff for /imach/src/imach.c between versions 1.237 and 1.250

version 1.237, 2016/08/26 09:20:19 version 1.250, 2016/09/08 16:07:27
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.250  2016/09/08 16:07:27  brouard
     Summary: continue
   
     Revision 1.249  2016/09/07 17:14:18  brouard
     Summary: Starting values from frequencies
   
     Revision 1.248  2016/09/07 14:10:18  brouard
     *** empty log message ***
   
     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    Revision 1.237  2016/08/26 09:20:19  brouard
   Summary: to valgrind    Summary: to valgrind
   
Line 1029  double dval; Line 1070  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 1145  double *Tvalsel; /**< Selected modality Line 1187  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 1154  int *Ndum; /** Freq of modality (tricode Line 1198  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 2057  void powell(double p[], double **xi, int Line 2100  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 2092  void powell(double p[], double **xi, int Line 2135  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 2170  void powell(double p[], double **xi, int Line 2246  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 2216  void powell(double p[], double **xi, int Line 2293  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 2496  Earliest age to start was %d-%d=%d, ncvl Line 2573  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 2519  Earliest age to start was %d-%d=%d, ncvl Line 2596  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 2555  Earliest age to start was %d-%d=%d, ncvl Line 2633  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 2610  Earliest age to start was %d-%d=%d, ncvl Line 2722  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 2743  double **pmij(double **ps, double *cov, Line 2860  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 3126  double func( double *x) Line 3244  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 3145  double func( double *x) Line 3263  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 3159  double func( double *x) Line 3278  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 3441  double funcone( double *x) Line 3563  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 3468  double funcone( double *x) Line 3591  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 3521  double funcone( double *x) Line 3649  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 3901  double hessij( double x[], double **hess Line 4029  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 4058  void pstamp(FILE *fichier) Line 4186  void pstamp(FILE *fichier)
 }  }
   
 /************ Frequencies ********************/  /************ Frequencies ********************/
 void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \  void  freqsummary(char fileres[], double p[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
                   int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[], \                    int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[], \
                   int firstpass,  int lastpass, int stepm, int weightopt, char model[])                    int firstpass,  int lastpass, int stepm, int weightopt, char model[])
 {  /* Some frequencies */  {  /* Some frequencies as well as proposing some starting values */
       
   int i, m, jk, j1, bool, z1,j, k, iv;    int i, m, jk, j1, bool, z1,j, k, iv, jj=0;
   int iind=0, iage=0;    int iind=0, iage=0;
   int mi; /* Effective wave */    int mi; /* Effective wave */
   int first;    int first;
Line 4090  void  freqsummary(char fileres[], int ia Line 4218  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 4100  void  freqsummary(char fileres[], int ia Line 4228  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 and dummy covariate value at beginning 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;    
   
   /* 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 
   */    */
     dateintsum=0;
     k2cpt=0;
   
     for (j = 0; j <= cptcoveff; j+=cptcoveff){   
     first=1;
   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 4158  Title=%s <br>Datafile=%s Firstpass=%d La Line 4288  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 */
     for (iind=1; iind<=imx; iind++) { /* For each individual iind */      for (iind=1; iind<=imx; iind++) { /* For each individual iind */
       bool=1;        bool=1;
         if(j !=0){
       if(anyvaryingduminmodel==0){ /* If All fixed covariates */        if(anyvaryingduminmodel==0){ /* If All fixed covariates */
         if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */          if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
           /* for (z1=1; z1<= nqfveff; z1++) {   */            /* for (z1=1; z1<= nqfveff; z1++) {   */
           /*   meanq[z1]+=coqvar[Tvar[z1]][iind];  /\* Computes mean of quantitative with selected filter *\/ */            /*   meanq[z1]+=coqvar[Tvar[z1]][iind];  /\* Computes mean of quantitative with selected filter *\/ */
           /* } */            /* } */
           for (z1=1; z1<=cptcoveff; z1++) {              for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */
             /* if(Tvaraff[z1] ==-20){ */              /* if(Tvaraff[z1] ==-20){ */
             /*   /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */              /*   /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
             /* }else  if(Tvaraff[z1] ==-10){ */              /* }else  if(Tvaraff[z1] ==-10){ */
             /*   /\* sumnew+=coqvar[z1][iind]; *\/ */              /*   /\* sumnew+=coqvar[z1][iind]; *\/ */
             /* }else  */              /* }else  */
             if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){              if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ /* for combination j1 of covariates */
               /* Tests if this individual iind responded to j1 (V4=1 V3=0) */                /* Tests if this individual iind responded to combination j1 (V4=1 V3=0) */
               bool=0;                bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */
               /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n",                 /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", 
                  bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),                   bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
                  j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/                   j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
Line 4186  Title=%s <br>Datafile=%s Firstpass=%d La Line 4318  Title=%s <br>Datafile=%s Firstpass=%d La
           } /* end z1 */            } /* end z1 */
         } /* cptcovn > 0 */          } /* cptcovn > 0 */
       } /* end any */        } /* end any */
         }/* end j==0 */
       if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */        if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */
         /* for(m=firstpass; m<=lastpass; m++){ */          /* for(m=firstpass; m<=lastpass; m++){ */
         for(mi=1; mi<wav[iind];mi++){ /* For that wave */          for(mi=1; mi<wav[iind];mi++){ /* For that wave */
           m=mw[mi][iind];            m=mw[mi][iind];
             if(j!=0){
           if(anyvaryingduminmodel==1){ /* Some are varying covariates */            if(anyvaryingduminmodel==1){ /* Some are varying covariates */
             for (z1=1; z1<=cptcoveff; z1++) {              for (z1=1; z1<=cptcoveff; z1++) {
               if( Fixed[Tmodelind[z1]]==1){                if( Fixed[Tmodelind[z1]]==1){
                 iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;                  iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
                 if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality */                  if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality. If covariate's 
                   bool=0;                                                                                    value is -1, we don't select. It differs from the 
                                                                                     constant and age model which counts them. */
                     bool=0; /* not selected */
               }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */                }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */
                 if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) {                  if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) {
                   bool=0;                    bool=0;
Line 4203  Title=%s <br>Datafile=%s Firstpass=%d La Line 4339  Title=%s <br>Datafile=%s Firstpass=%d La
               }                }
             }              }
           }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop  */            }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop  */
             } /* end j==0 */
           /* bool =0 we keep that guy which corresponds to the combination of dummy values */            /* bool =0 we keep that guy which corresponds to the combination of dummy values */
           if(bool==1){            if(bool==1){
             /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]              /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]
Line 4222  Title=%s <br>Datafile=%s Firstpass=%d La Line 4359  Title=%s <br>Datafile=%s Firstpass=%d La
                 if(s[m][iind]==-1)                  if(s[m][iind]==-1)
                   printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));                    printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));
                 freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */                  freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
                   /* if((int)agev[m][iind] == 55) */
                   /*   printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */
                 /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */                  /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */
                 freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */                  freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
               }                }
             } /* end if between passes */                } /* end if between passes */  
             if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) {              if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99) && (j==0)) {
               dateintsum=dateintsum+k2;                dateintsum=dateintsum+k2; /* on all covariates ?*/
               k2cpt++;                k2cpt++;
               /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */                /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
             }              }
           } /* end bool 2 */            }else{
               bool=1;
             }/* end bool 2 */
         } /* end m */          } /* end m */
       } /* end bool */        } /* end bool */
     } /* 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 && j!=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 4276  Title=%s <br>Datafile=%s Firstpass=%d La Line 4423  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 4401  Title=%s <br>Datafile=%s Firstpass=%d La Line 4548  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 */
     if(j==0){ /* We can estimate starting values from the occurences in each case */
       printf("#Freqsummary\n");
       fprintf(ficlog,"\n");
       for(i=1,jk=1; i <=nlstate; i++){
         for(k=1; k <=(nlstate+ndeath); k++){
           if (k != i) {
             printf("%d%d ",i,k);
             fprintf(ficlog,"%d%d ",i,k);
             for(jj=1; jj <=ncovmodel; jj++){
               if(jj==1){
                 printf("%12.7f ln(%12.1f/%12.1f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
                 fprintf(ficlog,"%12.7f ln(%12.1f/%12.1f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
               }
               /* printf("%12.7f )", param[i][jj][k]); */
               /* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */
               jk++; 
             }
             printf("\n");
             fprintf(ficlog,"\n");
           }
         }
       }
       printf("#Freqsummary\n");
       fprintf(ficlog,"\n");
       for(jk=-1; jk <=nlstate+ndeath; jk++){
         for(m=-1; m <=nlstate+ndeath; m++){
           /* param[i]|j][k]= freq[jk][m][iagemax+3] */
             printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]);
             fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]);
           /* if(freq[jk][m][iage] !=0 ) { /\* minimizing output *\/ */
           /*   printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */
           /*   fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */
           /* } */
         }
       } /* end loop jk */
       printf("\n");
       fprintf(ficlog,"\n");
     } /* if j=0 */
     } /* end j */
   dateintmean=dateintsum/k2cpt;     dateintmean=dateintsum/k2cpt; 
             
   fclose(ficresp);    fclose(ficresp);
   fclose(ficresphtm);    fclose(ficresphtm);
   fclose(ficresphtmfr);    fclose(ficresphtmfr);
Line 4737  void  concatwav(int wav[], int **dh, int Line 4923  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 4918  void  concatwav(int wav[], int **dh, int Line 5104  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 5283  void  concatwav(int wav[], int **dh, int Line 5469  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 5352  void  concatwav(int wav[], int **dh, int Line 5546  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 5384  void  concatwav(int wav[], int **dh, int Line 5578  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 5461  void  concatwav(int wav[], int **dh, int Line 5655  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 5556  void  concatwav(int wav[], int **dh, int Line 5750  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 5642  void  concatwav(int wav[], int **dh, int Line 5838  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 6046  void printinghtml(char fileresu[], char Line 6244  void printinghtml(char fileresu[], char
    jj1=0;     jj1=0;
   
    for(nres=1; nres <= nresult; nres++) /* For each resultline */     for(nres=1; nres <= nresult; nres++) /* For each resultline */
    for(k1=1; k1<=m;k1++){     for(k1=1; k1<=m;k1++){ /* For each combination of covariate */
      if(TKresult[nres]!= k1)       if(TKresult[nres]!= k1)
        continue;         continue;
   
Line 6074  void printinghtml(char fileresu[], char Line 6272  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 6177  See page 'Matrix of variance-covariance Line 6375  See page 'Matrix of variance-covariance
   
    jj1=0;     jj1=0;
   
    for(nres=1; nres <= nresult; nres++) /* For each resultline */     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)       if(TKresult[nres]!= k1)
        continue;         continue;
Line 6201  See page 'Matrix of variance-covariance Line 6399  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 6270  void printinggnuplot(char fileresu[], ch Line 6469  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)
             continue;
           /* We are interested in selected combination by the resultline */
           /* 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);
           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 */
             /* 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]; /* 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 */
             /* printf(" 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 */
             /* printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */
             fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
           }       
           /* printf("\n#\n"); */
           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-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
           fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);
           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);
         
           for (i=1; i<= nlstate ; i ++) {
             if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
             else        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);
           for (i=1; i<= nlstate ; i ++) {
             if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
             else 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); 
           for (i=1; i<= nlstate ; i ++) {
             if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
             else 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));
           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,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1, nres in 2 to be fixed */
             if(cptcoveff ==0){
               fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line lt 3",    2+(cpt-1),  cpt );
             }else{
               kl=0;
               for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */
                 lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
                 /* 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];
                 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' 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 */
     } /* cpt */
   
     
     /*2 eme*/
     for (k1=1; k1<= m ; k1 ++){  
       for(nres=1; nres <= nresult; nres++){ /* For each resultline */
       if(TKresult[nres]!= k1)        if(TKresult[nres]!= k1)
         continue;          continue;
       /* We are interested in selected combination by the resultline */        fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
       printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt);        for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
       fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files  and live state =%d ", cpt);          lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
       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 */  
         /* 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];
         /* 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);  
         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]);          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");  
       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,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);  
       fprintf(ficgp,"set xlabel \"Age\" \n\  
 set ylabel \"Probability\" \n             \  
 set ter svg size 640, 480\n                                             \  
 plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);  
         
       for (i=1; i<= nlstate ; i ++) {  
         if (i==cpt) fprintf(ficgp," %%lf (%%lf)");  
         else        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);  
       for (i=1; i<= nlstate ; i ++) {  
         if (i==cpt) fprintf(ficgp," %%lf (%%lf)");  
         else 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);   
       for (i=1; i<= nlstate ; i ++) {  
         if (i==cpt) fprintf(ficgp," %%lf (%%lf)");  
         else 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));  
       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,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */  
         if(cptcoveff ==0){  
           fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ",  2+(cpt-1),  cpt );  
         }else{  
           kl=0;  
           for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */  
             lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */  
             /* 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];  
             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++;  
             }  
           } /* end covariate */  
         } /* end if no covariate */  
       } /* end if backcast */  
       fprintf(ficgp,"\nset out \n");  
     } /* k1 */  
   } /* cpt */  
   
     
   /*2 eme*/  
   for (k1=1; k1<= m ; k1 ++)    
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */  
     if(TKresult[nres]!= k1)  
       continue;  
     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);        fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1,nres);
     for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/        for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
       if(vpopbased==0)          if(vpopbased==0)
         fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);            fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
       else          else
         fprintf(ficgp,"\nreplot ");            fprintf(ficgp,"\nreplot ");
       for (i=1; i<= nlstate+1 ; i ++) {          for (i=1; i<= nlstate+1 ; i ++) {
         k=2*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);            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 ++) {            for (j=1; j<= nlstate+1 ; j ++) {
           if (j==i) fprintf(ficgp," %%lf (%%lf)");              if (j==i) fprintf(ficgp," %%lf (%%lf)");
           else fprintf(ficgp," %%*lf (%%*lf)");              else fprintf(ficgp," %%*lf (%%*lf)");
         }               }   
         if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);            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);            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);            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 ++) {            for (j=1; j<= nlstate+1 ; j ++) {
           if (j==i) fprintf(ficgp," %%lf (%%lf)");              if (j==i) fprintf(ficgp," %%lf (%%lf)");
           else fprintf(ficgp," %%*lf (%%*lf)");              else fprintf(ficgp," %%*lf (%%*lf)");
         }               }   
         fprintf(ficgp,"\" t\"\" w l lt 0,");            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);            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 ++) {            for (j=1; j<= nlstate+1 ; j ++) {
           if (j==i) fprintf(ficgp," %%lf (%%lf)");              if (j==i) fprintf(ficgp," %%lf (%%lf)");
           else fprintf(ficgp," %%*lf (%%*lf)");              else fprintf(ficgp," %%*lf (%%*lf)");
         }               }   
         if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");            if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
         else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");            else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
       } /* state */          } /* state */
     } /* vpopbased */        } /* vpopbased */
     fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */        fprintf(ficgp,"\nset out;set out \"%s_%d-%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */
       } /* end nres */
   } /* k1 end 2 eme*/    } /* k1 end 2 eme*/
                   
                   
   /*3eme*/    /*3eme*/
   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]!= k)        if(TKresult[nres]!= k1)
       continue;  
   
     for (cpt=1; cpt<= nlstate ; cpt ++) {  
       fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files:  combination=%d state=%d",k1, cpt);  
       for (k=1; k<=cptcoveff; k++){    /* For each covariate dummy combination 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(k=1; k <= ncovds; k++){ */  
       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 ++) {
           fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files:  combination=%d state=%d",k1, cpt);
           for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
             lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
             /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
             /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
             /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
             vlv= nbcode[Tvaraff[k]][lv];
             fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
           }
           for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
             fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
           }       
           fprintf(ficgp,"\n#\n");
           if(invalidvarcomb[k1]){
             fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
             continue;
           }
                                                   
       /*       k=2+nlstate*(2*cpt-2); */          /*       k=2+nlstate*(2*cpt-2); */
       k=2+(nlstate+1)*(cpt-1);          k=2+(nlstate+1)*(cpt-1);
       fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);          fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres);
       fprintf(ficgp,"set ter svg size 640, 480\n\          fprintf(ficgp,"set ter svg size 640, 480\n\
 plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt);  plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt);
       /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);          /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
         for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");            for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
         fprintf(ficgp,"\" t \"e%d1\" w l",cpt);            fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
         fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);            fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
         for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");            for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
         fprintf(ficgp,"\" t \"e%d1\" w l",cpt);            fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
                                                                   
       */          */
       for (i=1; i< nlstate ; i ++) {          for (i=1; i< nlstate ; i ++) {
         fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);            fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
         /*      fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/            /*    fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
                                                                   
       }           } 
       fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);          fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
     }        }
   }      } /* end nres */
     } /* end kl 3eme */
       
   /* 4eme */    /* 4eme */
   /* Survival functions (period) from state i in state j by initial state i */    /* Survival functions (period) from state i in state j by initial state i */
   for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */    for (k1=1; k1<=m; k1++){    /* For each covariate and each value */
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */      for(nres=1; nres <= nresult; nres++){ /* For each resultline */
     if(TKresult[nres]!= k1)        if(TKresult[nres]!= k1)
       continue;  
   
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */  
       fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);  
       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 life state cpt*/
                                   fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);          for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\            lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
 set ter svg size 640, 480\n                                                                                                                                                                                     \            /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
 unset log y\n                                                                                                                                                                                                                                           \            /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
 plot [%.f:%.f]  ", ageminpar, agemaxpar);            /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
       k=3;            vlv= nbcode[Tvaraff[k]][lv];
       for (i=1; i<= nlstate ; i ++){            fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
         if(i==1){  
           fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));  
         }else{  
           fprintf(ficgp,", '' ");  
         }          }
         l=(nlstate+ndeath)*(i-1)+1;          for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
         fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);            fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
         for (j=2; j<= nlstate+ndeath ; j ++)          }       
           fprintf(ficgp,"+$%d",k+l+j-1);          fprintf(ficgp,"\n#\n");
         fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt);          if(invalidvarcomb[k1]){
       } /* nlstate */            fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
       fprintf(ficgp,"\nset out\n");            continue;
     } /* end cpt state*/           }
   } /* end covariate */          
                   fprintf(ficgp,"\nset out \"%s_%d-%d-%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 */
Line 6591  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6786  plot [%.f:%.f]  ", ageminpar, agemaxpar)
         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 6639  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6832  plot [%.f:%.f]  ", ageminpar, agemaxpar)
           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 6693  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6884  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 6764  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6953  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 6782  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6971  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 6797  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6987  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,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year \n");      fprintf(ficgp,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year\n");
     fprintf(ficgp,"#model=%s \n",model);      fprintf(ficgp,"#model=%s \n",model);
     fprintf(ficgp,"# ng=%d\n",ng);      fprintf(ficgp,"# Type of graphic ng=%d\n",ng);
     fprintf(ficgp,"#   jk=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */      fprintf(ficgp,"#   jk=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */
     for(jk=1; jk <=m; jk++)  /* For each combination of covariate */      for(jk=1; jk <=m; jk++)  /* For each combination of covariate */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */      for(nres=1; nres <= nresult; nres++){ /* For each resultline */
Line 6810  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 7000  plot [%.f:%.f]  ", ageminpar, agemaxpar)
         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");
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng);        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 6856  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 7046  plot [%.f:%.f]  ", ageminpar, agemaxpar)
               /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */                /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
               if(j==Tage[ij]) { /* Product by age */                if(j==Tage[ij]) { /* Product by age */
                 if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */                  if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
                   if(Dummy[j]==0){                    if(DummyV[j]==0){
                     fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);;                      fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);;
                   }else{ /* quantitative */                    }else{ /* quantitative */
                     fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */                      fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */
Line 6867  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 7057  plot [%.f:%.f]  ", ageminpar, agemaxpar)
               }else if(j==Tprod[ijp]) { /* */                 }else if(j==Tprod[ijp]) { /* */ 
                 /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */                  /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
                 if(ijp <=cptcovprod) { /* Product */                  if(ijp <=cptcovprod) { /* Product */
                   if(Dummy[Tvard[ijp][1]]==0){/* Vn is dummy */                    if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */
                     if(Dummy[Tvard[ijp][2]]==0){/* Vn and Vm are 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,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]]);                        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 */                      }else{ /* Vn is dummy and Vm is quanti */
Line 6876  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 7066  plot [%.f:%.f]  ", ageminpar, agemaxpar)
                       fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);                        fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
                     }                      }
                   }else{ /* Vn*Vm Vn is quanti */                    }else{ /* Vn*Vm Vn is quanti */
                     if(Dummy[Tvard[ijp][2]]==0){                      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]]);                        fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]);
                     }else{ /* Both quanti */                      }else{ /* Both quanti */
                       fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);                        fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
                     }                      }
                   }                    }
                     ijp++;
                 }                  }
               } else{  /* simple covariate */                } else{  /* simple covariate */
                 /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(jk,j)]); /\* Valgrind bug nbcode *\/ */                  /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(jk,j)]); /\* Valgrind bug nbcode *\/ */
Line 7175  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 7366  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 7785  int readdata(char datafile[], int firsto Line 7975  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 7822  int readdata(char datafile[], int firsto Line 8034  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 8207  int decodemodel( char model[], int lasto Line 8419  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 8413  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8625  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 8438  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8650  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 8452  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8664  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 8463  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8675  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 8479  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8691  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 8491  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8703  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 8693  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8905  Dummy[k] 0=dummy (0 1), 1 quantitative (
 }  }
   
 int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )  int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )
 {  {/* Check ages at death */
   int i, m;    int i, m;
   int firstone=0;    int firstone=0;
       
Line 9028  int prevalence_limit(double *p, double * Line 9240  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 9123  int back_prevalence_limit(double *p, dou Line 9337  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 9277  int hPijx(double *p, int bage, int fage) Line 9491  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 9305  int hPijx(double *p, int bage, int fage) Line 9519  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 9729  int main(int argc, char *argv[]) Line 9949  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 9737  int main(int argc, char *argv[]) Line 9963  int main(int argc, char *argv[])
     fclose (ficlog);      fclose (ficlog);
     goto end;      goto end;
     exit(0);      exit(0);
     } else if(mle==-2) { /* Guessing from means */
       prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
       printf(" You chose mle=-2, look at file %s for a template of covariance matrix \n",filereso);
       fprintf(ficlog," You chose mle=-2, look at file %s for a template of covariance matrix \n",filereso);
      
   }  else if(mle==-5) { /* Main Wizard */    }  else if(mle==-5) { /* Main Wizard */
     prwizard(ncovmodel, nlstate, ndeath, model, ficparo);      prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
     printf(" You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso);      printf(" You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
Line 10060  Please run with mle=-1 to get a correct Line 10291  Please run with mle=-1 to get a correct
   */    */
   
   concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);    concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
   /* */    /* Concatenates waves */
     
   free_vector(moisdc,1,n);    free_vector(moisdc,1,n);
   free_vector(andc,1,n);    free_vector(andc,1,n);
Line 10241  Title=%s <br>Datafile=%s Firstpass=%d La Line 10472  Title=%s <br>Datafile=%s Firstpass=%d La
   /* Calculates basic frequencies. Computes observed prevalence at single age     /* Calculates basic frequencies. Computes observed prevalence at single age 
                  and for any valid combination of covariates                   and for any valid combination of covariates
      and prints on file fileres'p'. */       and prints on file fileres'p'. */
   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \    freqsummary(fileres, p, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \
               firstpass, lastpass,  stepm,  weightopt, model);                firstpass, lastpass,  stepm,  weightopt, model);
   
   fprintf(fichtm,"\n");    fprintf(fichtm,"\n");
Line 10263  Interval (in months) between two waves: Line 10494  Interval (in months) between two waves:
   /* For mortality only */    /* For mortality only */
   if (mle==-3){    if (mle==-3){
     ximort=matrix(1,NDIM,1,NDIM);       ximort=matrix(1,NDIM,1,NDIM); 
                 for(i=1;i<=NDIM;i++)      for(i=1;i<=NDIM;i++)
                         for(j=1;j<=NDIM;j++)        for(j=1;j<=NDIM;j++)
                                 ximort[i][j]=0.;          ximort[i][j]=0.;
     /*     ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */      /*     ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */
     cens=ivector(1,n);      cens=ivector(1,n);
     ageexmed=vector(1,n);      ageexmed=vector(1,n);
Line 10501  Please run with mle=-1 to get a correct Line 10732  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 10751  Please run with mle=-1 to get a correct Line 10986  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 10767  Please run with mle=-1 to get a correct Line 11003  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 10785  Please run with mle=-1 to get a correct Line 11022  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 10803  Please run with mle=-1 to get a correct Line 11043  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 10818  Please run with mle=-1 to get a correct Line 11062  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 10953  Please run with mle=-1 to get a correct Line 11197  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 11034  Please run with mle=-1 to get a correct Line 11280  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 11065  Please run with mle=-1 to get a correct Line 11311  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 11242  Please run with mle=-1 to get a correct Line 11489  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.237  
changed lines
  Added in v.1.250

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