Diff for /imach/src/imach.c between versions 1.207 and 1.212

version 1.207, 2015/10/27 17:36:57 version 1.212, 2015/11/21 12:47:24
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.212  2015/11/21 12:47:24  brouard
     Summary: minor typo
   
     Revision 1.211  2015/11/21 12:41:11  brouard
     Summary: 0.98r3 with some graph of projected cross-sectional
   
     Author: Nicolas Brouard
   
     Revision 1.210  2015/11/18 17:41:20  brouard
     Summary: Start working on projected prevalences
   
     Revision 1.209  2015/11/17 22:12:03  brouard
     Summary: Adding ftolpl parameter
     Author: N Brouard
   
     We had difficulties to get smoothed confidence intervals. It was due
     to the period prevalence which wasn't computed accurately. The inner
     parameter ftolpl is now an outer parameter of the .imach parameter
     file after estepm. If ftolpl is small 1.e-4 and estepm too,
     computation are long.
   
     Revision 1.208  2015/11/17 14:31:57  brouard
     Summary: temporary
   
   Revision 1.207  2015/10/27 17:36:57  brouard    Revision 1.207  2015/10/27 17:36:57  brouard
   *** empty log message ***    *** empty log message ***
   
Line 741  typedef struct { Line 765  typedef struct {
 #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */  #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
 #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */  #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */
 #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1  #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
   /*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/
   #define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 
 #define MAXN 20000  #define MAXN 20000
 #define YEARM 12. /**< Number of months per year */  #define YEARM 12. /**< Number of months per year */
 #define AGESUP 130  #define AGESUP 130
Line 1985  double **prevalim(double **prlim, int nl Line 2011  double **prevalim(double **prlim, int nl
   /* 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;
   double min, max, maxmin, maxmax,sumnew=0.;    double *min, *max, *meandiff, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */    /* double **matprod2(); */ /* test */
   double **out, cov[NCOVMAX+1], **pmij();    double **out, cov[NCOVMAX+1], **pmij();
   double **newm;    double **newm;
   double agefin, delaymax=100 ; /* Max number of years to converge */    double agefin, delaymax=200. ; /* 100 Max number of years to converge */
   int ncvloop=0;    int ncvloop=0;
       
     min=vector(1,nlstate);
     max=vector(1,nlstate);
     meandiff=vector(1,nlstate);
   
   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);
Line 2029  double **prevalim(double **prlim, int nl Line 2059  double **prevalim(double **prlim, int nl
           
     savm=oldm;      savm=oldm;
     oldm=newm;      oldm=newm;
     maxmax=0.;  
     for(j=1;j<=nlstate;j++){      for(j=1; j<=nlstate; j++){
       min=1.;        max[j]=0.;
       max=0.;        min[j]=1.;
       for(i=1; i<=nlstate; i++) {      }
         sumnew=0;      for(i=1;i<=nlstate;i++){
         for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];        sumnew=0;
         for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];
         for(j=1; j<=nlstate; j++){ 
         prlim[i][j]= newm[i][j]/(1-sumnew);          prlim[i][j]= newm[i][j]/(1-sumnew);
         max=FMAX(max,prlim[i][j]);          max[j]=FMAX(max[j],prlim[i][j]);
         min=FMIN(min,prlim[i][j]);          min[j]=FMIN(min[j],prlim[i][j]);
         printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min);  
       }        }
       maxmin=(max-min)/(max+min)*2;      }
       maxmax=FMAX(maxmax,maxmin);  
       maxmax=0.;
       for(j=1; j<=nlstate; j++){
         meandiff[j]=(max[j]-min[j])/(max[j]+min[j])*2.; /* mean difference for each column */
         maxmax=FMAX(maxmax,meandiff[j]);
         /* printf(" age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, j, meandiff[j],(int)agefin, j, max[j], j, min[j],maxmax); */
     } /* j loop */      } /* j loop */
     *ncvyear= (int)age- (int)agefin;      *ncvyear= (int)age- (int)agefin;
     printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear);      /* printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */
     if(maxmax < ftolpl){      if(maxmax < ftolpl){
       /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */        /* printf("maxmax=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */
         free_vector(min,1,nlstate);
         free_vector(max,1,nlstate);
         free_vector(meandiff,1,nlstate);
       return prlim;        return prlim;
     }      }
   } /* age loop */    } /* age loop */
   printf("Warning: the stable prevalence at age %d did not converge with the required precision %g > ftolpl=%g. \n\      /* After some age loop it doesn't converge */
 Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);    printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
 /* 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); */  Earliest 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); */
     free_vector(min,1,nlstate);
     free_vector(max,1,nlstate);
     free_vector(meandiff,1,nlstate);
     
   return prlim; /* should not reach here */    return prlim; /* should not reach here */
 }  }
   
Line 2686  void likelione(FILE *ficres,double p[], Line 2730  void likelione(FILE *ficres,double p[],
       fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with optimized parameters mle = %d.",mle);        fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with optimized parameters mle = %d.",mle);
     fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));      fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
           
         
       for (k=1; k<= nlstate ; k++) {
         fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
   <img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k);
       }
     fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \      fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \
 <img src=\"%s-ori.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));  <img src=\"%s-ori.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
     fprintf(fichtm,"<br>- and by state of destination <a href=\"%s-dest.png\">%s-dest.png</a><br> \      fprintf(fichtm,"<br>- and by state of destination <a href=\"%s-dest.png\">%s-dest.png</a><br> \
 <img src=\"%s-dest.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));  <img src=\"%s-dest.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
     fflush(fichtm);      fflush(fichtm);
         
     for (k=1; k<= nlstate ; k++) {  
       fprintf(fichtm,"<br>- Probability p%dj by origin %d and destination j <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \  
 <img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k);  
     }  
   }    }
   return;    return;
 }  }
Line 2960  double hessij( double x[], double **hess Line 3004  double hessij( double x[], double **hess
   double p2[MAXPARM+1];    double p2[MAXPARM+1];
   int k, kmax=1;    int k, kmax=1;
   double v1, v2, cv12, lc1, lc2;    double v1, v2, cv12, lc1, lc2;
   
     int firstime=0;
       
   fx=func(x);    fx=func(x);
   for (k=1; k<=kmax; k=k+10) {    for (k=1; k<=kmax; k=k+10) {
Line 2981  double hessij( double x[], double **hess Line 3027  double hessij( double x[], double **hess
     k4=func(p2)-fx;      k4=func(p2)-fx;
     res=(k1-k2-k3+k4)/4.0/delti[thetai]/k/delti[thetaj]/k/2.; /* Because of L not 2*L */      res=(k1-k2-k3+k4)/4.0/delti[thetai]/k/delti[thetaj]/k/2.; /* Because of L not 2*L */
     if(k1*k2*k3*k4 <0.){      if(k1*k2*k3*k4 <0.){
         firstime=1;
       kmax=kmax+10;        kmax=kmax+10;
       if(kmax >=10){      }
       if(kmax >=10 || firstime ==1){
       printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol);        printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol);
       fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol);        fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; 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);
       }  
     }      }
 #ifdef DEBUGHESSIJ  #ifdef DEBUGHESSIJ
     v1=hess[thetai][thetai];      v1=hess[thetai][thetai];
Line 3187  void  freqsummary(char fileres[], int ia Line 3234  void  freqsummary(char fileres[], int ia
       k2cpt=0;        k2cpt=0;
       for (i=1; i<=imx; i++) {        for (i=1; i<=imx; i++) {
         bool=1;          bool=1;
         if  (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */          if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
           for (z1=1; z1<=cptcoveff; z1++)                   for (z1=1; z1<=cptcoveff; z1++)       
             if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){              if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
                 /* Tests if the value of each of the covariates of i is equal to filter j1 */                  /* Tests if the value of each of the covariates of i is equal to filter j1 */
Line 3197  void  freqsummary(char fileres[], int ia Line 3244  void  freqsummary(char fileres[], int ia
                 j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/                  j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
               /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/                /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
             }               } 
         }          } /* cptcovn > 0 */
     
         if (bool==1){          if (bool==1){
           for(m=firstpass; m<=lastpass; m++){            for(m=firstpass; m<=lastpass; m++){
Line 3211  void  freqsummary(char fileres[], int ia Line 3258  void  freqsummary(char fileres[], int ia
                 freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i];                  freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i];
               }                }
                               
               if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3))) {                if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3)) && (anint[m][i]!=9999) && (mint[m][i]!=99)) {
                 dateintsum=dateintsum+k2;                  dateintsum=dateintsum+k2;
                 k2cpt++;                  k2cpt++;
                   /* printf("i=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",i, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
               }                }
               /*}*/                /*}*/
           }            } /* end m */
         }          } /* end bool */
       } /* end i */        } /* end i = 1 to imx */
                 
       /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/        /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
       pstamp(ficresp);        pstamp(ficresp);
Line 3304  void  freqsummary(char fileres[], int ia Line 3352  void  freqsummary(char fileres[], int ia
         if(first==1)          if(first==1)
           printf("Others in log...\n");            printf("Others in log...\n");
         fprintf(ficlog,"\n");          fprintf(ficlog,"\n");
       }        } /* end loop i */
       /*}*/        /*}*/
   }    } /* end j1 */
   dateintmean=dateintsum/k2cpt;     dateintmean=dateintsum/k2cpt; 
     
   fclose(ficresp);    fclose(ficresp);
Line 3993  void cvevsij(double ***eij, double x[], Line 4041  void cvevsij(double ***eij, double x[],
 }  }
   
 /************ Variance ******************/  /************ Variance ******************/
  void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[])   void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[])
 {  {
   /* Variance of health expectancies */    /* Variance of health expectancies */
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/    /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
Line 4047  void cvevsij(double ***eij, double x[], Line 4095  void cvevsij(double ***eij, double x[],
     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);      fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
   }    }
   printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);    printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
    
   fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);    fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
   pstamp(ficresprobmorprev);    pstamp(ficresprobmorprev);
   fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);    fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
Line 4058  void cvevsij(double ***eij, double x[], Line 4105  void cvevsij(double ***eij, double x[],
       fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j);        fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j);
   }      }  
   fprintf(ficresprobmorprev,"\n");    fprintf(ficresprobmorprev,"\n");
     
   fprintf(ficgp,"\n# Routine varevsij");    fprintf(ficgp,"\n# Routine varevsij");
   fprintf(ficgp,"\nunset title \n");    fprintf(ficgp,"\nunset title \n");
 /* fprintf(fichtm, "#Local time at start: %s", strstart);*/  /* fprintf(fichtm, "#Local time at start: %s", strstart);*/
Line 4095  void cvevsij(double ***eij, double x[], Line 4143  void cvevsij(double ***eij, double x[],
   /* For example we decided to compute the life expectancy with the smallest unit */    /* For example we decided to compute the life expectancy with the smallest unit */
   /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.     /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. 
      nhstepm is the number of hstepm from age to agelim        nhstepm is the number of hstepm from age to agelim 
      nstepm is the number of stepm from age to agelin.        nstepm is the number of stepm from age to agelim. 
      Look at function hpijx to understand why (it is linked to memory size questions) */       Look at function hpijx to understand why because of memory size limitations, 
   /* We decided (b) to get a life expectancy respecting the most precise curvature of the       we decided (b) to get a life expectancy respecting the most precise curvature of the
      survival function given by stepm (the optimization length). Unfortunately it       survival function given by stepm (the optimization length). Unfortunately it
      means that if the survival funtion is printed every two years of age and if       means that if the survival funtion is printed every two years of age and if
      you sum them up and add 1 year (area under the trapezoids) you won't get the same        you sum them up and add 1 year (area under the trapezoids) you won't get the same 
Line 4118  void cvevsij(double ***eij, double x[], Line 4166  void cvevsij(double ***eij, double x[],
       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);
       }        }
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);    
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
   
       if (popbased==1) {        if (popbased==1) {
         if(mobilav ==0){          if(mobilav ==0){
Line 4131  void cvevsij(double ***eij, double x[], Line 4179  void cvevsij(double ***eij, double x[],
         }          }
       }        }
       
         hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  /* Returns p3mat[i][j][h] for h=1 to nhstepm */
       for(j=1; j<= nlstate; j++){        for(j=1; j<= nlstate; j++){
         for(h=0; h<=nhstepm; h++){          for(h=0; h<=nhstepm; h++){
           for(i=1, gp[h][j]=0.;i<=nlstate;i++)            for(i=1, gp[h][j]=0.;i<=nlstate;i++)
             gp[h][j] += prlim[i][i]*p3mat[i][j][h];              gp[h][j] += prlim[i][i]*p3mat[i][j][h];
         }          }
       }        }
       /* This for computing probability of death (h=1 means        /* Next for computing probability of death (h=1 means
          computed over hstepm matrices product = hstepm*stepm months)            computed over hstepm matrices product = hstepm*stepm months) 
          as a weighted average of prlim.           as a weighted average of prlim.
       */        */
Line 4149  void cvevsij(double ***eij, double x[], Line 4198  void cvevsij(double ***eij, double x[],
   
       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);
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);    
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear, ij);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij);
     
       if (popbased==1) {        if (popbased==1) {
         if(mobilav ==0){          if(mobilav ==0){
Line 4162  void cvevsij(double ***eij, double x[], Line 4211  void cvevsij(double ***eij, double x[],
         }          }
       }        }
   
         hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
   
       for(j=1; j<= nlstate; j++){  /* Sum of wi * eij = e.j */        for(j=1; j<= nlstate; j++){  /* Sum of wi * eij = e.j */
         for(h=0; h<=nhstepm; h++){          for(h=0; h<=nhstepm; h++){
           for(i=1, gm[h][j]=0.;i<=nlstate;i++)            for(i=1, gm[h][j]=0.;i<=nlstate;i++)
Line 4224  void cvevsij(double ***eij, double x[], Line 4275  void cvevsij(double ***eij, double x[],
         varppt[j][i]=doldmp[j][i];          varppt[j][i]=doldmp[j][i];
     /* end ppptj */      /* end ppptj */
     /*  x centered again */      /*  x centered again */
     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);    
     prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyear,ij);      prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij);
     
     if (popbased==1) {      if (popbased==1) {
       if(mobilav ==0){        if(mobilav ==0){
Line 4241  void cvevsij(double ***eij, double x[], Line 4292  void cvevsij(double ***eij, double x[],
        computed over hstepm (estepm) matrices product = hstepm*stepm months)          computed over hstepm (estepm) matrices product = hstepm*stepm months) 
        as a weighted average of prlim.         as a weighted average of prlim.
     */      */
       hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);  
     for(j=nlstate+1;j<=nlstate+ndeath;j++){      for(j=nlstate+1;j<=nlstate+ndeath;j++){
       for(i=1,gmp[j]=0.;i<= nlstate; i++)         for(i=1,gmp[j]=0.;i<= nlstate; i++) 
         gmp[j] += prlim[i][i]*p3mat[i][j][1];           gmp[j] += prlim[i][i]*p3mat[i][j][1]; 
Line 4303  void cvevsij(double ***eij, double x[], Line 4355  void cvevsij(double ***eij, double x[],
 }  /* end varevsij */  }  /* end varevsij */
   
 /************ Variance of prevlim ******************/  /************ Variance of prevlim ******************/
  void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, char strstart[])   void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[])
 {  {
   /* Variance of prevalence limit  for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/    /* Variance of prevalence limit  for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/    /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/
Line 4313  void cvevsij(double ***eij, double x[], Line 4365  void cvevsij(double ***eij, double x[],
   double *xp;    double *xp;
   double *gp, *gm;    double *gp, *gm;
   double **gradg, **trgradg;    double **gradg, **trgradg;
     double **mgm, **mgp;
   double age,agelim;    double age,agelim;
   int theta;    int theta;
       
Line 4335  void cvevsij(double ***eij, double x[], Line 4388  void cvevsij(double ***eij, double x[],
     if (stepm >= YEARM) hstepm=1;      if (stepm >= YEARM) hstepm=1;
     nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */      nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
     gradg=matrix(1,npar,1,nlstate);      gradg=matrix(1,npar,1,nlstate);
       mgp=matrix(1,npar,1,nlstate);
       mgm=matrix(1,npar,1,nlstate);
     gp=vector(1,nlstate);      gp=vector(1,nlstate);
     gm=vector(1,nlstate);      gm=vector(1,nlstate);
   
Line 4342  void cvevsij(double ***eij, double x[], Line 4397  void cvevsij(double ***eij, double x[],
       for(i=1; i<=npar; i++){ /* Computes gradient */        for(i=1; i<=npar; i++){ /* Computes gradient */
         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,ncvyear,ij);        if((int)age==79 ||(int)age== 80 ||(int)age== 81 )
       for(i=1;i<=nlstate;i++)          prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
         else
           prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
         for(i=1;i<=nlstate;i++){
         gp[i] = prlim[i][i];          gp[i] = prlim[i][i];
               mgp[theta][i] = prlim[i][i];
         }
       for(i=1; i<=npar; i++) /* Computes gradient */        for(i=1; i<=npar; i++) /* Computes gradient */
         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,ncvyear,ij);        if((int)age==79 ||(int)age== 80 ||(int)age== 81 )
       for(i=1;i<=nlstate;i++)          prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
         else
           prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
         for(i=1;i<=nlstate;i++){
         gm[i] = prlim[i][i];          gm[i] = prlim[i][i];
           mgm[theta][i] = prlim[i][i];
         }
       for(i=1;i<=nlstate;i++)        for(i=1;i<=nlstate;i++)
         gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta];          gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta];
         /* gradg[theta][2]= -gradg[theta][1]; */ /* For testing if nlstate=2 */
     } /* End theta */      } /* End theta */
   
     trgradg =matrix(1,nlstate,1,npar);      trgradg =matrix(1,nlstate,1,npar);
Line 4361  void cvevsij(double ***eij, double x[], Line 4425  void cvevsij(double ***eij, double x[],
     for(j=1; j<=nlstate;j++)      for(j=1; j<=nlstate;j++)
       for(theta=1; theta <=npar; theta++)        for(theta=1; theta <=npar; theta++)
         trgradg[j][theta]=gradg[theta][j];          trgradg[j][theta]=gradg[theta][j];
       /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */
       /*   printf("\nmgm mgp %d ",(int)age); */
       /*   for(j=1; j<=nlstate;j++){ */
       /*  printf(" %d ",j); */
       /*  for(theta=1; theta <=npar; theta++) */
       /*    printf(" %d %lf %lf",theta,mgm[theta][j],mgp[theta][j]); */
       /*  printf("\n "); */
       /*   } */
       /* } */
       /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */
       /*   printf("\n gradg %d ",(int)age); */
       /*   for(j=1; j<=nlstate;j++){ */
       /*  printf("%d ",j); */
       /*  for(theta=1; theta <=npar; theta++) */
       /*    printf("%d %lf ",theta,gradg[theta][j]); */
       /*  printf("\n "); */
       /*   } */
       /* } */
   
     for(i=1;i<=nlstate;i++)      for(i=1;i<=nlstate;i++)
       varpl[i][(int)age] =0.;        varpl[i][(int)age] =0.;
     if((int)age==67 ||(int)age== 66 ){      if((int)age==79 ||(int)age== 80  ||(int)age== 81){
     matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);      matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);
     matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);      matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);
     }else{      }else{
Line 4380  void cvevsij(double ***eij, double x[], Line 4462  void cvevsij(double ***eij, double x[],
     fprintf(ficresvpl,"\n");      fprintf(ficresvpl,"\n");
     free_vector(gp,1,nlstate);      free_vector(gp,1,nlstate);
     free_vector(gm,1,nlstate);      free_vector(gm,1,nlstate);
       free_matrix(mgm,1,npar,1,nlstate);
       free_matrix(mgp,1,npar,1,nlstate);
     free_matrix(gradg,1,npar,1,nlstate);      free_matrix(gradg,1,npar,1,nlstate);
     free_matrix(trgradg,1,nlstate,1,npar);      free_matrix(trgradg,1,nlstate,1,npar);
   } /* End age */    } /* End age */
Line 4730  To be simple, these graphs help to under Line 4814  To be simple, these graphs help to under
 void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \  void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \
                   int lastpass, int stepm, int weightopt, char model[],\                    int lastpass, int stepm, int weightopt, char model[],\
                   int imx,int jmin, int jmax, double jmeanint,char rfileres[],\                    int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
                   int popforecast, int estepm ,\                    int popforecast, int prevfcast, int estepm ,          \
                   double jprev1, double mprev1,double anprev1, \                    double jprev1, double mprev1,double anprev1, \
                   double jprev2, double mprev2,double anprev2){                    double jprev2, double mprev2,double anprev2){
   int jj1, k1, i1, cpt;    int jj1, k1, i1, cpt;
Line 4748  void printinghtml(char fileresu[], char Line 4832  void printinghtml(char fileresu[], char
  - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",   - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
            subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_"));             subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_"));
    fprintf(fichtm,"\     fprintf(fichtm,"\
  - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \   - (a) Life expectancies by health status at initial age, e<sub>i.</sub> (b) health expectancies by health status at initial age, e<sub>ij</sub> . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
    <a href=\"%s\">%s</a> <br>\n",     <a href=\"%s\">%s</a> <br>\n",
            estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_"));             estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_"));
    fprintf(fichtm,"\     if(prevfcast==1){
  - Population projections by age and states: \       fprintf(fichtm,"\
    - Prevalence projections by age and states:                            \
    <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_"));     <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_"));
      }
   
 fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");  fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
   
Line 4773  fprintf(fichtm," \n<ul><li><b>Graphs</b> Line 4859  fprintf(fichtm," \n<ul><li><b>Graphs</b>
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }       }
      /* aij, bij */       /* aij, bij */
      fprintf(fichtm,"<br>- Logit model, 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.svg\">%s_%d-1.svg</a><br> \
 <img src=\"%s_%d-1.svg\">",subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);  <img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
      /* Pij */       /* Pij */
      fprintf(fichtm,"<br>\n- Pij 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.svg\">%s_%d-2.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.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);     
      /* Quasi-incidences */       /* Quasi-incidences */
      fprintf(fichtm,"<br>\n- Iij 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 hPij \   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: hPij/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.svg\">%s_%d-3.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.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); 
      /* Survival functions (period) in state j */       /* Survival functions (period) in state j */
      for(cpt=1; cpt<=nlstate;cpt++){       for(cpt=1; cpt<=nlstate;cpt++){
Line 4791  divided by h: hPij/h : <a href=\"%s_%d-3 Line 4877  divided by h: hPij/h : <a href=\"%s_%d-3
      }       }
      /* 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 any different live states 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.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);
      }       }
Line 4800  divided by h: hPij/h : <a href=\"%s_%d-3 Line 4886  divided by h: hPij/h : <a href=\"%s_%d-3
        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.svg\">%s%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.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);
      }       }
       if(prevfcast==1){
         /* Projection of prevalence up to period (stable) prevalence in each health state */
         for(cpt=1; cpt<=nlstate;cpt++){
           fprintf(fichtm,"<br>\n- Projection of prevalece 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> \
   <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1);
         }
       }
   
      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.svg\">%s_%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.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1);
Line 4887  true period expectancies (those weighted Line 4981  true period expectancies (those weighted
 }  }
   
 /******************* Gnuplot file **************/  /******************* Gnuplot file **************/
 void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){      void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, char pathc[], double p[]){
   
   char dirfileres[132],optfileres[132];    char dirfileres[132],optfileres[132];
   int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;    int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
     int lv=0, vlv=0, kl=0;
   int ng=0;    int ng=0;
   int vpopbased;    int vpopbased;
 /*   if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */  /*   if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
Line 4935  void printinggnuplot(char fileresu[], ch Line 5030  void printinggnuplot(char fileresu[], ch
   strcpy(dirfileres,optionfilefiname);    strcpy(dirfileres,optionfilefiname);
   strcpy(optfileres,"vpl");    strcpy(optfileres,"vpl");
  /* 1eme*/   /* 1eme*/
   fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files\n");    for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */
   for (cpt=1; cpt<= nlstate ; cpt ++) {      for (k1=1; k1<= m ; k1 ++) { /* For each combination of covariate */
     for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */        /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
         fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files ");
         for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
           lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
           /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
           /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
           /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
           vlv= nbcode[Tvaraff[lv]][lv];
           fprintf(ficgp," V%d=%d ",k,vlv);
         }
         fprintf(ficgp,"\n#\n");
   
      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
      fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);       fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
      fprintf(ficgp,"set xlabel \"Age\" \n\       fprintf(ficgp,"set xlabel \"Age\" \n\
Line 4964  plot [%.f:%.f] \"%s\" every :::%d::%d u Line 5070  plot [%.f:%.f] \"%s\" every :::%d::%d u
     } /* k1 */      } /* k1 */
   } /* cpt */    } /* cpt */
   /*2 eme*/    /*2 eme*/
   fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n");  
   for (k1=1; k1<= m ; k1 ++) {     for (k1=1; k1<= m ; k1 ++) { 
         fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
         for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
           lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
           /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
           /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
           /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
           vlv= nbcode[Tvaraff[lv]][lv];
           fprintf(ficgp," V%d=%d ",k,vlv);
         }
         fprintf(ficgp,"\n#\n");
   
     fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);      fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
     for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/      for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
       if(vpopbased==0)        if(vpopbased==0)
Line 4998  plot [%.f:%.f] \"%s\" every :::%d::%d u Line 5114  plot [%.f:%.f] \"%s\" every :::%d::%d u
     } /* 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.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
   } /* k1 */    } /* k1 */
   
   
   /*3eme*/    /*3eme*/
     
   for (k1=1; k1<= m ; k1 ++) {     for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<= nlstate ; cpt ++) {      for (cpt=1; cpt<= nlstate ; cpt ++) {
         fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files:  cov=%d state=%d",k1, cpt);
         for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
           lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
           /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
           /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
           /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
           vlv= nbcode[Tvaraff[lv]][lv];
           fprintf(ficgp," V%d=%d ",k,vlv);
         }
         fprintf(ficgp,"\n#\n");
   
       /*       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.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
Line 5027  plot [%.f:%.f] \"%s\" every :::%d::%d u Line 5155  plot [%.f:%.f] \"%s\" every :::%d::%d u
   /* Survival functions (period) from state i in state j by initial state i */    /* Survival functions (period) from state i in state j by initial state i */
   for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */    for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
       k=3;        fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
       fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'lij' files, cov=%d state=%d",k1, cpt);        for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
           lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
           /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
           /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
           /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
           vlv= nbcode[Tvaraff[lv]][lv];
           fprintf(ficgp," V%d=%d ",k,vlv);
         }
         fprintf(ficgp,"\n#\n");
   
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);        fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
 set ter svg size 640, 480\n\  set ter svg size 640, 480\n\
 unset log y\n\  unset log y\n\
 plot [%.f:%.f]  ", ageminpar, agemaxpar);  plot [%.f:%.f]  ", ageminpar, agemaxpar);
         k=3;
       for (i=1; i<= nlstate ; i ++){        for (i=1; i<= nlstate ; i ++){
         if(i==1)          if(i==1)
           fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));            fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
Line 5052  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 5190  plot [%.f:%.f]  ", ageminpar, agemaxpar)
   /* Survival functions (period) from state i in state j by final state j */    /* Survival functions (period) from state i in state j by final state j */
   for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */    for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
       k=3;  
       fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);        fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
         for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
           lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
           /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
           /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
           /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
           vlv= nbcode[Tvaraff[lv]][lv];
           fprintf(ficgp," V%d=%d ",k,vlv);
         }
         fprintf(ficgp,"\n#\n");
   
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);        fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
 set ter svg size 640, 480\n\  set ter svg size 640, 480\n\
 unset log y\n\  unset log y\n\
 plot [%.f:%.f]  ", ageminpar, agemaxpar);  plot [%.f:%.f]  ", ageminpar, agemaxpar);
         k=3;
       for (j=1; j<= nlstate ; j ++){ /* Lived in state j */        for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
         if(j==1)          if(j==1)
           fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));            fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
Line 5084  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 5232  plot [%.f:%.f]  ", ageminpar, agemaxpar)
   } /* end covariate */      } /* end covariate */  
   
   /* CV preval stable (period) for each covariate */    /* CV preval stable (period) for each covariate */
   for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */    for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
       k=3;        fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
       fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt);        for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
           lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
           /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
           /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
           /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
           vlv= nbcode[Tvaraff[lv]][lv];
           fprintf(ficgp," V%d=%d ",k,vlv);
         }
         fprintf(ficgp,"\n#\n");
   
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);        fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
 set ter svg size 640, 480\n\  set ter svg size 640, 480\n\
 unset log y\n\  unset log y\n\
 plot [%.f:%.f]  ", ageminpar, agemaxpar);  plot [%.f:%.f]  ", ageminpar, agemaxpar);
         k=3; /* Offset */
       for (i=1; i<= nlstate ; i ++){        for (i=1; i<= nlstate ; i ++){
         if(i==1)          if(i==1)
           fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));            fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
Line 5108  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 5266  plot [%.f:%.f]  ", ageminpar, agemaxpar)
     } /* end cpt state*/       } /* end cpt state*/ 
   } /* end covariate */      } /* end covariate */  
   
     if(prevfcast==1){
     /* Projection from cross-sectional to stable (period) for each covariate */
   
       for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
         for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
           fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
           for (k=1; k<=cptcoveff; k++){    /* For each correspondig covariate value  */
             lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
             /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
             /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
             /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
             vlv= nbcode[Tvaraff[lv]][lv];
             fprintf(ficgp," V%d=%d ",k,vlv);
           }
           fprintf(ficgp,"\n#\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,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
   set ter svg size 640, 480\n\
   unset log y\n\
   plot [%.f:%.f]  ", ageminpar, agemaxpar);
           for (i=1; i<= nlstate+1 ; i ++){  /* nlstate +1 p11 p21 p.1 */
             /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
             /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */   
             /*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
             /*#   1       2   3    4    5      6  7   8   9   10   11 12  13   14  15 */   
             if(i==1){
               fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_"));
             }else{
               fprintf(ficgp,",\\\n '' ");
             }
             if(cptcoveff ==0){ /* No covariate */
               fprintf(ficgp," u 2:("); /* Age is in 2 */
               /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
               /*#   1       2   3   4   5  6    7  8   9   10  11  12  13  14  15  16  17  18 */
               if(i==nlstate+1)
                 fprintf(ficgp," $%d/(1.-$%d)) t 'p.%d' with line ", \
                           2+(cpt-1)*(nlstate+1)+1+(i-1),  2+1+(i-1)+(nlstate+1)*nlstate,cpt );
               else
                 fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \
                         2+(cpt-1)*(nlstate+1)+1+(i-1),  2+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
             }else{
               fprintf(ficgp,"u 6:(("); /* Age is in 6 */
               /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
               /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */   
               kl=0;
               for (k=1; k<=cptcoveff; k++){    /* For each covariate  */
                 lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
                 /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
                 /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
                 /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
                 vlv= nbcode[Tvaraff[lv]][lv];
                 kl++;
                 /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
                 /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ 
                 /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ 
                 /* ''  u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
                 if(k==cptcoveff)
                   if(i==nlstate+1)
                     fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \
                             6+(cpt-1)*(nlstate+1)+1+(i-1),  6+1+(i-1)+(nlstate+1)*nlstate,cpt );
                   else
                     fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \
                             6+(cpt-1)*(nlstate+1)+1+(i-1),  6+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
                 else{
                   fprintf(ficgp,"$%d==%d && $%d==%d && ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv]);
                   kl++;
                 }
               } /* end covariate */
             } /* end if covariate */
           } /* nlstate */
           fprintf(ficgp,"\nset out\n");
         } /* end cpt state*/
       } /* end covariate */
     } /* End if prevfcast */
   
   
   /* proba elementaires */    /* proba elementaires */
   fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");    fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
   for(i=1,jk=1; i <=nlstate; i++){    for(i=1,jk=1; i <=nlstate; i++){
Line 5299  void prevforecast(char fileres[], double Line 5535  void prevforecast(char fileres[], double
   char fileresf[FILENAMELENGTH];    char fileresf[FILENAMELENGTH];
   
   agelim=AGESUP;    agelim=AGESUP;
     /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
        in each health status at the date of interview (if between dateprev1 and dateprev2).
        We still use firstpass and lastpass as another selection.
     */
   prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);    prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
     
   strcpy(fileresf,"F_");     strcpy(fileresf,"F_"); 
Line 5349  void prevforecast(char fileres[], double Line 5589  void prevforecast(char fileres[], double
   for(cptcov=1, k=0;cptcov<=i1;cptcov++){    for(cptcov=1, k=0;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){      for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
       k=k+1;        k=k+1;
       fprintf(ficresf,"\n#******");        fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#");
       for(j=1;j<=cptcoveff;j++) {        for(j=1;j<=cptcoveff;j++) {
         fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);          fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       }        }
       fprintf(ficresf,"******\n");        fprintf(ficresf," yearproj age");
       fprintf(ficresf,"# Covariate valuofcovar yearproj age");  
       for(j=1; j<=nlstate+ndeath;j++){         for(j=1; j<=nlstate+ndeath;j++){ 
         for(i=1; i<=nlstate;i++)                        for(i=1; i<=nlstate;i++)              
           fprintf(ficresf," p%d%d",i,j);            fprintf(ficresf," p%d%d",i,j);
Line 6510  void syscompilerinfo(int logged) Line 6749  void syscompilerinfo(int logged)
   
  }   }
   
  int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyear){   int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
   /*--------------- Prevalence limit  (period or stable prevalence) --------------*/    /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
   int i, j, k, i1 ;    int i, j, k, i1 ;
   /* double ftolpl = 1.e-10; */    /* double ftolpl = 1.e-10; */
Line 6566  void syscompilerinfo(int logged) Line 6805  void syscompilerinfo(int logged)
                   
         for (age=agebase; age<=agelim; age++){          for (age=agebase; age<=agelim; age++){
         /* for (age=agebase; age<=agebase; age++){ */          /* for (age=agebase; age<=agebase; age++){ */
           prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k);            prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
           fprintf(ficrespl,"%.0f ",age );            fprintf(ficrespl,"%.0f ",age );
           for(j=1;j<=cptcoveff;j++)            for(j=1;j<=cptcoveff;j++)
             fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);              fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
Line 6575  void syscompilerinfo(int logged) Line 6814  void syscompilerinfo(int logged)
             tot +=  prlim[i][i];              tot +=  prlim[i][i];
             fprintf(ficrespl," %.5f", prlim[i][i]);              fprintf(ficrespl," %.5f", prlim[i][i]);
           }            }
           fprintf(ficrespl," %.3f %d\n", tot, *ncvyear);            fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
         } /* Age */          } /* Age */
         /* was end of cptcod */          /* was end of cptcod */
     } /* cptcov */      } /* cptcov */
Line 6668  int main(int argc, char *argv[]) Line 6907  int main(int argc, char *argv[])
 #endif  #endif
   int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);    int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
   int i,j, k, n=MAXN,iter=0,m,size=100, cptcod;    int i,j, k, n=MAXN,iter=0,m,size=100, cptcod;
   int ncvyearnp=0;    int ncvyear=0; /* Number of years needed for the period prevalence to converge */
   int *ncvyear=&ncvyearnp; /* Number of years needed for the period prevalence to converge */  
   int jj, ll, li, lj, lk;    int jj, ll, li, lj, lk;
   int numlinepar=0; /* Current linenumber of parameter file */    int numlinepar=0; /* Current linenumber of parameter file */
   int num_filled;    int num_filled;
Line 6929  int main(int argc, char *argv[]) Line 7167  int main(int argc, char *argv[])
   if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \    if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \
                         &ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){                          &ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){
     if (num_filled != 8) {      if (num_filled != 8) {
       printf("Not 8\n");        printf("Not 8 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n");
         printf("but line=%s\n",line);
     }      }
     printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt);      printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt);
   }    }
   /* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */    /* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */
   ftolpl=6.e-3; /* 6.e-3 make convergences in less than 80 loops for the prevalence limit */    /*ftolpl=6.e-4; *//* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
   /* Third parameter line */    /* Third parameter line */
   while(fgets(line, MAXLINE, ficpar)) {    while(fgets(line, MAXLINE, ficpar)) {
     /* If line starts with a # it is a comment */      /* If line starts with a # it is a comment */
Line 7303  Please run with mle=-1 to get a correct Line 7542  Please run with mle=-1 to get a correct
   Ndum =ivector(-1,NCOVMAX);      Ndum =ivector(-1,NCOVMAX);  
   if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */    if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */
     tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */      tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */
   /* Nbcode gives the value of the lth modality of jth covariate, in    /* Nbcode gives the value of the lth modality (currently 1 to 2) of jth covariate, in
      V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/       V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/
   /* 1 to ncodemax[j] is the maximum value of this jth covariate */    /* 1 to ncodemax[j] which is the maximum value of this jth covariate */
   
   /*  codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */    /*  codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
   /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/    /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/
   /* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/    /* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/
     /* nbcode[Tvaraff[j]][codtabm(h,j)]) : if there are only 2 modalities for a covariate j, 
      * codtabm(h,j) gives its value classified at position h and nbcode gives how it is coded 
      * (currently 0 or 1) in the data.
      * In a loop on h=1 to 2**k, and a loop on j (=1 to k), we get the value of 
      * corresponding modality (h,j).
      */
   
   h=0;    h=0;
   
   
Line 7319  Please run with mle=-1 to get a correct Line 7565  Please run with mle=-1 to get a correct
   m=pow(2,cptcoveff);    m=pow(2,cptcoveff);
     
           /**< codtab(h,k)  k   = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1            /**< codtab(h,k)  k   = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1
            * For k=4 covariates, h goes from 1 to 2**k             * For k=4 covariates, h goes from 1 to m=2**k
            * codtabm(h,k)=  1 & (h-1) >> (k-1) ;             * codtabm(h,k)=  (1 & (h-1) >> (k-1)) + 1;
              * #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
            *     h\k   1     2     3     4             *     h\k   1     2     3     4
            *______________________________               *______________________________  
            *     1 i=1 1 i=1 1 i=1 1 i=1 1             *     1 i=1 1 i=1 1 i=1 1 i=1 1
Line 7340  Please run with mle=-1 to get a correct Line 7587  Please run with mle=-1 to get a correct
            *    15 i=8 1     2     2     2             *    15 i=8 1     2     2     2
            *    16     2     2     2     2             *    16     2     2     2     2
            */             */
     /* How to do the opposite? From combination h (=1 to 2**k) how to get the value on the covariates? */
        /* from h=5 and m, we get then number of covariates k=log(m)/log(2)=4
        * and the value of each covariate?
        * V1=1, V2=1, V3=2, V4=1 ?
        * h-1=4 and 4 is 0100 or reverse 0010, and +1 is 1121 ok.
        * h=6, 6-1=5, 5 is 0101, 1010, 2121, V1=2nd, V2=1st, V3=2nd, V4=1st.
        * In order to get the real value in the data, we use nbcode
        * nbcode[Tvar[3][2nd]]=1 and nbcode[Tvar[4][1]]=0
        * We are keeping this crazy system in order to be able (in the future?) 
        * to have more than 2 values (0 or 1) for a covariate.
        * #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
        * h=6, k=2? h-1=5=0101, reverse 1010, +1=2121, k=2nd position: value is 1: codtabm(6,2)=1
        *              bbbbbbbb
        *              76543210     
        *   h-1        00000101 (6-1=5)
        *(h-1)>>(k-1)= 00000001 >> (2-1) = 1 right shift
        *           &
        *     1        00000001 (1)
        *              00000001        = 1 & ((h-1) >> (k-1))
        *          +1= 00000010 =2 
        *
        * h=14, k=3 => h'=h-1=13, k'=k-1=2
        *          h'      1101 =2^3+2^2+0x2^1+2^0
        *    >>k'            11
        *          &   00000001
        *            = 00000001
        *      +1    = 00000010=2    =  codtabm(14,3)   
        * Reverse h=6 and m=16?
        * cptcoveff=log(16)/log(2)=4 covariate: 6-1=5=0101 reversed=1010 +1=2121 =>V1=2, V2=1, V3=2, V4=1.
        * for (j=1 to cptcoveff) Vj=decodtabm(j,h,cptcoveff)
        * decodtabm(h,j,cptcoveff)= (((h-1) >> (j-1)) & 1) +1 
        * decodtabm(h,j,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (j-1)) & 1) +1 : -1)
        * V3=decodtabm(14,3,2**4)=2
        *          h'=13   1101 =2^3+2^2+0x2^1+2^0
        *(h-1) >> (j-1)    0011 =13 >> 2
        *          &1 000000001
        *           = 000000001
        *         +1= 000000010 =2
        *                  2211
        *                  V1=1+1, V2=0+1, V3=1+1, V4=1+1
        *                  V3=2
        */
   
   /* /\* for(h=1; h <=100 ;h++){  *\/ */    /* /\* for(h=1; h <=100 ;h++){  *\/ */
   /*   /\* printf("h=%2d ", h); *\/ */    /*   /\* printf("h=%2d ", h); *\/ */
   /*    /\* for(k=1; k <=10; k++){ *\/ */    /*    /\* for(k=1; k <=10; k++){ *\/ */
Line 7870  Please run with mle=-1 to get a correct Line 8160  Please run with mle=-1 to get a correct
           
     fflush(ficlog);      fflush(ficlog);
     fflush(ficres);      fflush(ficres);
             while(fgets(line, MAXLINE, ficpar)) {
     while((c=getc(ficpar))=='#' && c!= EOF){      /* If line starts with a # it is a comment */
       ungetc(c,ficpar);      if (line[0] == '#') {
       fgets(line, MAXLINE, ficpar);        numlinepar++;
       fputs(line,stdout);        fputs(line,stdout);
       fputs(line,ficparo);        fputs(line,ficparo);
     }        fputs(line,ficlog);
     ungetc(c,ficpar);        continue;
       }else
         break;
     }
   
       /* while((c=getc(ficpar))=='#' && c!= EOF){ */
       /*   ungetc(c,ficpar); */
       /*   fgets(line, MAXLINE, ficpar); */
       /*   fputs(line,stdout); */
       /*   fputs(line,ficparo); */
       /* } */
       /* ungetc(c,ficpar); */
           
     estepm=0;      estepm=0;
     fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);      if((num_filled=sscanf(line,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm, &ftolpl)) !=EOF){
   
       if (num_filled != 6) {
         printf("Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n");
         printf("but line=%s\n",line);
         goto end;
       }
       printf("agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",ageminpar,agemaxpar, bage, fage, estepm, ftolpl);
     }
     /* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */
     /*ftolpl=6.e-4;*/ /* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
   
       /* fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); */
     if (estepm==0 || estepm < stepm) estepm=stepm;      if (estepm==0 || estepm < stepm) estepm=stepm;
     if (fage <= 2) {      if (fage <= 2) {
       bage = ageminpar;        bage = ageminpar;
Line 7888  Please run with mle=-1 to get a correct Line 8201  Please run with mle=-1 to get a correct
     }      }
           
     fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");      fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
     fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);      fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
     fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);      fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d, ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
   
     /* Other stuffs, more or less useful */          /* Other stuffs, more or less useful */    
     while((c=getc(ficpar))=='#' && c!= EOF){      while((c=getc(ficpar))=='#' && c!= EOF){
Line 7952  Please run with mle=-1 to get a correct Line 8265  Please run with mle=-1 to get a correct
 This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\  This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
 Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);  Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
     }else      }else
       printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);        printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, pathc,p);
           
     printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\      printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\
                  model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\                   model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,estepm, \
                  jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);                   jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
               
    /*------------ free_vector  -------------*/     /*------------ free_vector  -------------*/
Line 7979  Please run with mle=-1 to get a correct Line 8292  Please run with mle=-1 to get a correct
     /*--------------- Prevalence limit  (period or stable prevalence) --------------*/      /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
     /*#include "prevlim.h"*/  /* Use ficrespl, ficlog */      /*#include "prevlim.h"*/  /* Use ficrespl, ficlog */
     prlim=matrix(1,nlstate,1,nlstate);      prlim=matrix(1,nlstate,1,nlstate);
     prevalence_limit(p, prlim,  ageminpar, agemaxpar, ftolpl, ncvyear);      prevalence_limit(p, prlim,  ageminpar, agemaxpar, ftolpl, &ncvyear);
     fclose(ficrespl);      fclose(ficrespl);
   
 #ifdef FREEEXIT2  #ifdef FREEEXIT2
Line 8042  Please run with mle=-1 to get a correct Line 8355  Please run with mle=-1 to get a correct
       printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);        printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
       fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0);        fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
     }      }
     printf("Computing Health Expectancies: result on file '%s' \n", filerese);      printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout);
     fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);      fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog);
     /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){      /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/        for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
                       
Line 8062  Please run with mle=-1 to get a correct Line 8375  Please run with mle=-1 to get a correct
       /*}*/        /*}*/
     }      }
     fclose(ficreseij);      fclose(ficreseij);
       printf("done evsij\n");fflush(stdout);
       fprintf(ficlog,"done evsij\n");fflush(ficlog);
   
     /*---------- Health expectancies and variances ------------*/      /*---------- Health expectancies and variances ------------*/
   
Line 8073  Please run with mle=-1 to get a correct Line 8387  Please run with mle=-1 to get a correct
       printf("Problem with total LE resultfile: %s\n", filerest);goto end;        printf("Problem with total LE resultfile: %s\n", filerest);goto end;
       fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;        fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
     }      }
     printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);       printf("Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(stdout);
     fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);       fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(ficlog);
   
   
     strcpy(fileresstde,"STDE_");      strcpy(fileresstde,"STDE_");
Line 8083  Please run with mle=-1 to get a correct Line 8397  Please run with mle=-1 to get a correct
       printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);        printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
       fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);        fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
     }      }
     printf("Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);      printf("  Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
     fprintf(ficlog,"Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);      fprintf(ficlog,"  Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
   
     strcpy(filerescve,"CVE_");      strcpy(filerescve,"CVE_");
     strcat(filerescve,fileresu);      strcat(filerescve,fileresu);
Line 8092  Please run with mle=-1 to get a correct Line 8406  Please run with mle=-1 to get a correct
       printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);        printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
       fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);        fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
     }      }
     printf("Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);      printf("    Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
     fprintf(ficlog,"Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);      fprintf(ficlog,"    Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
   
     strcpy(fileresv,"V_");      strcpy(fileresv,"V_");
     strcat(fileresv,fileresu);      strcat(fileresv,fileresu);
Line 8101  Please run with mle=-1 to get a correct Line 8415  Please run with mle=-1 to get a correct
       printf("Problem with variance resultfile: %s\n", fileresv);exit(0);        printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
       fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0);        fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0);
     }      }
     printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);      printf("      Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);fflush(stdout);
     fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);      fprintf(ficlog,"      Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);fflush(ficlog);
   
     /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){      /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/        for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
                       
     for (k=1; k <= (int) pow(2,cptcoveff); k++){      for (k=1; k <= (int) pow(2,cptcoveff); k++){
         fprintf(ficrest,"\n#****** ");        fprintf(ficrest,"\n#****** ");
         for(j=1;j<=cptcoveff;j++)         for(j=1;j<=cptcoveff;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)]);
         fprintf(ficrest,"******\n");        fprintf(ficrest,"******\n");
         
         fprintf(ficresstdeij,"\n#****** ");        fprintf(ficresstdeij,"\n#****** ");
         fprintf(ficrescveij,"\n#****** ");        fprintf(ficrescveij,"\n#****** ");
         for(j=1;j<=cptcoveff;j++) {        for(j=1;j<=cptcoveff;j++) {
           fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);          fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
           fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);          fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
         }        }
         fprintf(ficresstdeij,"******\n");        fprintf(ficresstdeij,"******\n");
         fprintf(ficrescveij,"******\n");        fprintf(ficrescveij,"******\n");
         
         fprintf(ficresvij,"\n#****** ");        fprintf(ficresvij,"\n#****** ");
         for(j=1;j<=cptcoveff;j++)         for(j=1;j<=cptcoveff;j++) 
           fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);          fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
         fprintf(ficresvij,"******\n");        fprintf(ficresvij,"******\n");
         
         eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);        eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
         oldm=oldms;savm=savms;        oldm=oldms;savm=savms;
         cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);          printf(" cvevsij %d, ",k);
         /*        fprintf(ficlog, " cvevsij %d, ",k);
          */        cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
         /* goto endfree; */        printf(" end cvevsij \n ");
          fprintf(ficlog, " end cvevsij \n ");
         vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);        
         pstamp(ficrest);        /*
          */
         /* goto endfree; */
         for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/        
           oldm=oldms;savm=savms; /* ZZ Segmentation fault */        vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
           cptcod= 0; /* To be deleted */        pstamp(ficrest);
           varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */        
           fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");        
           if(vpopbased==1)        for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
             fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);          oldm=oldms;savm=savms; /* ZZ Segmentation fault */
           else          cptcod= 0; /* To be deleted */
             fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");          printf("varevsij %d \n",vpopbased);
           fprintf(ficrest,"# Age popbased mobilav e.. (std) ");          fprintf(ficlog, "varevsij %d \n",vpopbased);
           for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);          varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
           fprintf(ficrest,"\n");          fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");
           /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */          if(vpopbased==1)
           epj=vector(1,nlstate+1);            fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
           for(age=bage; age <=fage ;age++){          else
             prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k); /*ZZ Is it the correct prevalim */            fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
             if (vpopbased==1) {          fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
               if(mobilav ==0){          for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
                 for(i=1; i<=nlstate;i++)          fprintf(ficrest,"\n");
                   prlim[i][i]=probs[(int)age][i][k];          /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
               }else{ /* mobilav */           epj=vector(1,nlstate+1);
                 for(i=1; i<=nlstate;i++)          printf("Computing age specific period (stable) prevalences in each health state \n");
                   prlim[i][i]=mobaverage[(int)age][i][k];          fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
               }          for(age=bage; age <=fage ;age++){
             }            prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */
                     if (vpopbased==1) {
             fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);              if(mobilav ==0){
             /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */                for(i=1; i<=nlstate;i++)
             /* printf(" age %4.0f ",age); */                  prlim[i][i]=probs[(int)age][i][k];
             for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){              }else{ /* mobilav */ 
               for(i=1, epj[j]=0.;i <=nlstate;i++) {                for(i=1; i<=nlstate;i++)
                 epj[j] += prlim[i][i]*eij[i][j][(int)age];                  prlim[i][i]=mobaverage[(int)age][i][k];
                 /*ZZZ  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/  
                 /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */  
               }  
               epj[nlstate+1] +=epj[j];  
             }              }
             /* printf(" age %4.0f \n",age); */            }
             
             for(i=1, vepp=0.;i <=nlstate;i++)            fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
               for(j=1;j <=nlstate;j++)            /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
                 vepp += vareij[i][j][(int)age];            /* printf(" age %4.0f ",age); */
             fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));            for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
             for(j=1;j <=nlstate;j++){              for(i=1, epj[j]=0.;i <=nlstate;i++) {
               fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));                epj[j] += prlim[i][i]*eij[i][j][(int)age];
                 /*ZZZ  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
                 /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
             }              }
             fprintf(ficrest,"\n");              epj[nlstate+1] +=epj[j];
             }
             /* printf(" age %4.0f \n",age); */
             
             for(i=1, vepp=0.;i <=nlstate;i++)
               for(j=1;j <=nlstate;j++)
                 vepp += vareij[i][j][(int)age];
             fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
             for(j=1;j <=nlstate;j++){
               fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
           }            }
             fprintf(ficrest,"\n");
         }          }
         free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);        } /* End vpopbased */
         free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);        free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
         free_vector(epj,1,nlstate+1);        free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
         free_vector(epj,1,nlstate+1);
         printf("done \n");fflush(stdout);
         fprintf(ficlog,"done\n");fflush(ficlog);
         
       /*}*/        /*}*/
     }      } /* End k */
     free_vector(weight,1,n);      free_vector(weight,1,n);
     free_imatrix(Tvard,1,NCOVMAX,1,2);      free_imatrix(Tvard,1,NCOVMAX,1,2);
     free_imatrix(s,1,maxwav+1,1,n);      free_imatrix(s,1,maxwav+1,1,n);
Line 8203  Please run with mle=-1 to get a correct Line 8529  Please run with mle=-1 to get a correct
     fclose(ficrescveij);      fclose(ficrescveij);
     fclose(ficresvij);      fclose(ficresvij);
     fclose(ficrest);      fclose(ficrest);
       printf("done Health expectancies\n");fflush(stdout);
       fprintf(ficlog,"done Health expectancies\n");fflush(ficlog);
     fclose(ficpar);      fclose(ficpar);
       
     /*------- Variance of period (stable) prevalence------*/         /*------- Variance of period (stable) prevalence------*/   
Line 8213  Please run with mle=-1 to get a correct Line 8541  Please run with mle=-1 to get a correct
       printf("Problem with variance of period (stable) prevalence  resultfile: %s\n", fileresvpl);        printf("Problem with variance of period (stable) prevalence  resultfile: %s\n", fileresvpl);
       exit(0);        exit(0);
     }      }
     printf("Computing Variance-covariance of period (stable) prevalence: file '%s' \n", fileresvpl);      printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout);
       fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog);
   
     /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){      /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/        for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
Line 8226  Please run with mle=-1 to get a correct Line 8555  Please run with mle=-1 to get a correct
               
         varpl=matrix(1,nlstate,(int) bage, (int) fage);          varpl=matrix(1,nlstate,(int) bage, (int) fage);
         oldm=oldms;savm=savms;          oldm=oldms;savm=savms;
         varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, k, strstart);          varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart);
         free_matrix(varpl,1,nlstate,(int) bage, (int)fage);          free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
       /*}*/        /*}*/
     }      }
   
     fclose(ficresvpl);      fclose(ficresvpl);
       printf("done variance-covariance of period prevalence\n");fflush(stdout);
       fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);
   
     /*---------- End : free ----------------*/      /*---------- End : free ----------------*/
     if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);      if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);

Removed from v.1.207  
changed lines
  Added in v.1.212


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