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

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

version 1.240, 2016/08/29 07:53:18 version 1.248, 2016/09/07 14:10:18
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     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    Revision 1.240  2016/08/29 07:53:18  brouard
   Summary: Better    Summary: Better
   
Line 2070  void powell(double p[], double **xi, int Line 2094  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,jk,k;    int i,ibig,j,jk,k; 
   double del,t,*pt,*ptt,*xit;    double del,t,*pt,*ptt,*xit;
Line 2111  void powell(double p[], double **xi, int Line 2135  void powell(double p[], double **xi, int
     printf("\n#model=  1      +     age ");      printf("\n#model=  1      +     age ");
     fprintf(ficlog,"\n#model=  1      +     age ");      fprintf(ficlog,"\n#model=  1      +     age ");
     if(nagesqr==1){      if(nagesqr==1){
         printf("  + age*age  ",Tvar[j]);          printf("  + age*age  ");
         fprintf(ficlog,"  + age*age  ",Tvar[j]);          fprintf(ficlog,"  + age*age  ");
     }      }
     for(j=1;j <=ncovmodel-2;j++){      for(j=1;j <=ncovmodel-2;j++){
       if(Typevar[j]==0) {        if(Typevar[j]==0) {
Line 2145  void powell(double p[], double **xi, int Line 2169  void powell(double p[], double **xi, int
         }          }
       }        }
     }      }
     if(*iter <=3){      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 2216  void powell(double p[], double **xi, int Line 2240  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 2542  Earliest age to start was %d-%d=%d, ncvl Line 2567  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 2565  Earliest age to start was %d-%d=%d, ncvl Line 2590  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 2601  Earliest age to start was %d-%d=%d, ncvl Line 2627  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 2656  Earliest age to start was %d-%d=%d, ncvl Line 2716  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 2789  double **pmij(double **ps, double *cov, Line 2854  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 3172  double func( double *x) Line 3238  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 3191  double func( double *x) Line 3257  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 3205  double func( double *x) Line 3272  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 3487  double funcone( double *x) Line 3557  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 3514  double funcone( double *x) Line 3585  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 3567  double funcone( double *x) Line 3643  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 3947  double hessij( double x[], double **hess Line 4023  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 4788  void  concatwav(int wav[], int **dh, int Line 4864  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 5411  void  concatwav(int wav[], int **dh, int Line 5487  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 5443  void  concatwav(int wav[], int **dh, int Line 5519  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 5520  void  concatwav(int wav[], int **dh, int Line 5596  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 5615  void  concatwav(int wav[], int **dh, int Line 5691  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 5701  void  concatwav(int wav[], int **dh, int Line 5779  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 6105  void printinghtml(char fileresu[], char Line 6185  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 6133  void printinghtml(char fileresu[], char Line 6213  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 6236  See page 'Matrix of variance-covariance Line 6316  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 6260  See page 'Matrix of variance-covariance Line 6340  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 6336  void printinggnuplot(char fileresu[], ch Line 6417  void printinggnuplot(char fileresu[], ch
         if(TKresult[nres]!= k1)          if(TKresult[nres]!= k1)
           continue;            continue;
         /* We are interested in selected combination by the resultline */          /* We are interested in selected combination by the resultline */
         printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt);          /* printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); */
         fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files  and live state =%d ", cpt);          fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files  and live state =%d ", cpt);
         for (k=1; k<=cptcoveff; k++){    /* For each covariate k get corresponding value lv for combination k1 */          for (k=1; k<=cptcoveff; k++){    /* For each covariate k get corresponding value lv for combination k1 */
           lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */            lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
Line 6345  void printinggnuplot(char fileresu[], ch Line 6426  void printinggnuplot(char fileresu[], ch
           /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */            /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
           vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */            vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */
           /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */            /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
           printf(" V%d=%d ",Tvaraff[k],vlv);            /* printf(" V%d=%d ",Tvaraff[k],vlv); */
           fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);            fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
         }          }
         for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */          for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
           printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);            /* printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */
           fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);            fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
         }                 }       
         printf("\n#\n");          /* printf("\n#\n"); */
         fprintf(ficgp,"\n#\n");          fprintf(ficgp,"\n#\n");
         if(invalidvarcomb[k1]){          if(invalidvarcomb[k1]){
           fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);             fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
           continue;            continue;
         }          }
               
         fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);          fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
         fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);          fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);
         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);          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 ++) {          for (i=1; i<= nlstate ; i ++) {
           if (i==cpt) fprintf(ficgp," %%lf (%%lf)");            if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
           else        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);          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 ++) {          for (i=1; i<= nlstate ; i ++) {
           if (i==cpt) fprintf(ficgp," %%lf (%%lf)");            if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
           else 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);           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 ++) {          for (i=1; i<= nlstate ; i ++) {
           if (i==cpt) fprintf(ficgp," %%lf (%%lf)");            if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
           else fprintf(ficgp," %%*lf (%%*lf)");            else fprintf(ficgp," %%*lf (%%*lf)");
Line 6380  void printinggnuplot(char fileresu[], ch Line 6461  void printinggnuplot(char fileresu[], ch
         fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));          fprintf(ficgp,"\" 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 */          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\" 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 */            fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1, nres in 2 to be fixed */
           if(cptcoveff ==0){            if(cptcoveff ==0){
             fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ",        2+(cpt-1),  cpt );              fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line lt 3",    2+(cpt-1),  cpt );
           }else{            }else{
             kl=0;              kl=0;
             for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */              for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */
Line 6397  void printinggnuplot(char fileresu[], ch Line 6478  void printinggnuplot(char fileresu[], ch
               /*6+1+(i-1)+(nlstate+1)*nlstate; 6+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*/                /* ''  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){                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], \                  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], \
                         4+(cpt-1),  cpt );  /* 4 or 6 ?*/                          2+cptcoveff*2+(cpt-1),  cpt );  /* 4 or 6 ?*/
               }else{                }else{
                 fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);                  fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
                 kl++;                  kl++;
Line 6437  void printinggnuplot(char fileresu[], ch Line 6518  void printinggnuplot(char fileresu[], ch
         continue;          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);
Line 6467  void printinggnuplot(char fileresu[], ch Line 6548  void printinggnuplot(char fileresu[], ch
           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 */      } /* end nres */
   } /* k1 end 2 eme*/    } /* k1 end 2 eme*/
                   
Line 6499  void printinggnuplot(char fileresu[], ch Line 6580  void printinggnuplot(char fileresu[], ch
                                                   
         /*       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);
Line 6545  plot [%.f:%.f] \"%s\" every :::%d::%d u Line 6626  plot [%.f:%.f] \"%s\" every :::%d::%d u
           continue;            continue;
         }          }
               
         fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);          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\          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);  set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
         k=3;          k=3;
Line 6591  set ter svg size 640, 480\nunset log y\n Line 6672  set ter svg size 640, 480\nunset log y\n
           continue;            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\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);  set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
         k=3;          k=3;
Line 6646  set ter svg size 640, 480\nunset log y\n Line 6727  set ter svg size 640, 480\nunset log y\n
         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\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);  set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
       k=3; /* Offset */        k=3; /* Offset */
Line 6692  set ter svg size 640, 480\nunset log y\n Line 6773  set ter svg size 640, 480\nunset log y\n
           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\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);  set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
         k=3; /* Offset */          k=3; /* Offset */
Line 6744  set ter svg size 640, 480\nunset log y\n Line 6825  set ter svg size 640, 480\nunset log y\n
         }          }
                   
         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\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);  set ter svg size 640, 480\nunset log y\nplot [%.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 */
Line 6860  set ter svg size 640, 480\nunset log y\n Line 6941  set ter svg size 640, 480\nunset log y\n
         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 8524  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8605  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 8535  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8616  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 8551  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8632  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 8765  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 8846  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 9238  int back_prevalence_limit(double *p, dou Line 9319  int back_prevalence_limit(double *p, dou
         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++)
Line 9809  int main(int argc, char *argv[]) Line 9890  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 9817  int main(int argc, char *argv[]) Line 9904  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 10140  Please run with mle=-1 to get a correct Line 10232  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 10343  Interval (in months) between two waves: Line 10435  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 10581  Please run with mle=-1 to get a correct Line 10673  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 11125  Please run with mle=-1 to get a correct Line 11221  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)]);
Removed from v.1.240  
changed lines
  Added in v.1.248

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