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

Diff for /imach/src/imach.c between versions 1.116 and 1.123

version 1.116, 2006/03/06 10:29:27 version 1.123, 2006/03/20 10:52:43
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.123  2006/03/20 10:52:43  brouard
     * imach.c (Module): <title> changed, corresponds to .htm file
     name. <head> headers where missing.
   
     * imach.c (Module): Weights can have a decimal point as for
     English (a comma might work with a correct LC_NUMERIC environment,
     otherwise the weight is truncated).
     Modification of warning when the covariates values are not 0 or
     1.
     Version 0.98g
   
     Revision 1.122  2006/03/20 09:45:41  brouard
     (Module): Weights can have a decimal point as for
     English (a comma might work with a correct LC_NUMERIC environment,
     otherwise the weight is truncated).
     Modification of warning when the covariates values are not 0 or
     1.
     Version 0.98g
   
     Revision 1.121  2006/03/16 17:45:01  lievre
     * imach.c (Module): Comments concerning covariates added
   
     * imach.c (Module): refinements in the computation of lli if
     status=-2 in order to have more reliable computation if stepm is
     not 1 month. Version 0.98f
   
     Revision 1.120  2006/03/16 15:10:38  lievre
     (Module): refinements in the computation of lli if
     status=-2 in order to have more reliable computation if stepm is
     not 1 month. Version 0.98f
   
     Revision 1.119  2006/03/15 17:42:26  brouard
     (Module): Bug if status = -2, the loglikelihood was
     computed as likelihood omitting the logarithm. Version O.98e
   
     Revision 1.118  2006/03/14 18:20:07  brouard
     (Module): varevsij Comments added explaining the second
     table of variances if popbased=1 .
     (Module): Covariances of eij, ekl added, graphs fixed, new html link.
     (Module): Function pstamp added
     (Module): Version 0.98d
   
     Revision 1.117  2006/03/14 17:16:22  brouard
     (Module): varevsij Comments added explaining the second
     table of variances if popbased=1 .
     (Module): Covariances of eij, ekl added, graphs fixed, new html link.
     (Module): Function pstamp added
     (Module): Version 0.98d
   
   Revision 1.116  2006/03/06 10:29:27  brouard    Revision 1.116  2006/03/06 10:29:27  brouard
   (Module): Variance-covariance wrong links and    (Module): Variance-covariance wrong links and
   varian-covariance of ej. is needed (Saito).    varian-covariance of ej. is needed (Saito).
Line 210 Line 259
   hPijx.    hPijx.
   
   Also this programme outputs the covariance matrix of the parameters but also    Also this programme outputs the covariance matrix of the parameters but also
   of the life expectancies. It also computes the stable prevalence.     of the life expectancies. It also computes the period (stable) prevalence. 
       
   Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).    Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).
            Institut national d'études démographiques, Paris.             Institut national d'études démographiques, Paris.
Line 239 Line 288
       begin-prev-date,...        begin-prev-date,...
   open gnuplot file    open gnuplot file
   open html file    open html file
   stable prevalence    period (stable) prevalence
    for age prevalim()     for age prevalim()
   h Pij x    h Pij x
   variance of p varprob    variance of p varprob
Line 251 Line 300
   varevsij()     varevsij() 
   if popbased==1 varevsij(,popbased)    if popbased==1 varevsij(,popbased)
   total life expectancies    total life expectancies
   Variance of stable prevalence    Variance of period (stable) prevalence
  end   end
 */  */
   
Line 311  extern int errno; Line 360  extern int errno;
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
   
 char version[]="Imach version 0.98c, February 2006, INED-EUROREVES ";  char version[]="Imach version 0.98g, March 2006, INED-EUROREVES-Institut de longevite ";
 char fullversion[]="$Revision$ $Date$";   char fullversion[]="$Revision$ $Date$"; 
   char strstart[80];
   char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */  int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
 int nvar;  int nvar;
 int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;  int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
Line 336  int **bh; /* bh[mi][i] is the bias (+ or Line 387  int **bh; /* bh[mi][i] is the bias (+ or
 double jmean; /* Mean space between 2 waves */  double jmean; /* Mean space between 2 waves */
 double **oldm, **newm, **savm; /* Working pointers to matrices */  double **oldm, **newm, **savm; /* Working pointers to matrices */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */  double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
 FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;  FILE *fic,*ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
 FILE *ficlog, *ficrespow;  FILE *ficlog, *ficrespow;
 int globpr; /* Global variable for printing or not */  int globpr; /* Global variable for printing or not */
 double fretone; /* Only one call to likelihood */  double fretone; /* Only one call to likelihood */
Line 350  FILE *ficresprobmorprev; Line 401  FILE *ficresprobmorprev;
 FILE *fichtm, *fichtmcov; /* Html File */  FILE *fichtm, *fichtmcov; /* Html File */
 FILE *ficreseij;  FILE *ficreseij;
 char filerese[FILENAMELENGTH];  char filerese[FILENAMELENGTH];
   FILE *ficresstdeij;
   char fileresstde[FILENAMELENGTH];
   FILE *ficrescveij;
   char filerescve[FILENAMELENGTH];
 FILE  *ficresvij;  FILE  *ficresvij;
 char fileresv[FILENAMELENGTH];  char fileresv[FILENAMELENGTH];
 FILE  *ficresvpl;  FILE  *ficresvpl;
 char fileresvpl[FILENAMELENGTH];  char fileresvpl[FILENAMELENGTH];
 char title[MAXLINE];  char title[MAXLINE];
 char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];  char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH];  char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH];
 char tmpout[FILENAMELENGTH],  tmpout2[FILENAMELENGTH];   char tmpout[FILENAMELENGTH],  tmpout2[FILENAMELENGTH]; 
 char command[FILENAMELENGTH];  char command[FILENAMELENGTH];
 int  outcmd=0;  int  outcmd=0;
Line 380  char strcurr[80], strfor[80]; Line 435  char strcurr[80], strfor[80];
   
 char *endptr;  char *endptr;
 long lval;  long lval;
   double dval;
   
 #define NR_END 1  #define NR_END 1
 #define FREE_ARG char*  #define FREE_ARG char*
Line 1069  void powell(double p[], double **xi, int Line 1125  void powell(double p[], double **xi, int
   }     } 
 }   } 
   
 /**** Prevalence limit (stable prevalence)  ****************/  /**** Prevalence limit (stable or period prevalence)  ****************/
   
 double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij)  double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij)
 {  {
Line 1362  double func( double *x) Line 1418  double func( double *x)
   
         } else if  (s2==-2) {          } else if  (s2==-2) {
           for (j=1,survp=0. ; j<=nlstate; j++)             for (j=1,survp=0. ; j<=nlstate; j++) 
             survp += out[s1][j];              survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
           lli= survp;            /*survp += out[s1][j]; */
         }            lli= log(survp);
           
         else if  (s2==-4) {  
           for (j=3,survp=0. ; j<=nlstate; j++)   
             survp += out[s1][j];  
           lli= survp;  
         }          }
                   
         else if  (s2==-5) {          else if  (s2==-4) { 
           for (j=1,survp=0. ; j<=2; j++)             for (j=3,survp=0. ; j<=nlstate; j++)  
             survp += out[s1][j];              survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
           lli= survp;            lli= log(survp); 
         }          } 
   
   
           else if  (s2==-5) { 
             for (j=1,survp=0. ; j<=2; j++)  
               survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
             lli= log(survp); 
           } 
           
         else{          else{
           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= (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 */
Line 1573  double funcone( double *x) Line 1629  double funcone( double *x)
        */         */
       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 (mle==1){        } else if  (s2==-2) {
           for (j=1,survp=0. ; j<=nlstate; j++) 
             survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
           lli= log(survp);
         }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 */
Line 1589  double funcone( double *x) Line 1649  double funcone( double *x)
       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,"%9d %6d %1d %1d %1d %1d %3d %10.6f %6.4f\          fprintf(ficresilk,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\
  %10.6f %10.6f %10.6f ", \   %11.6f %11.6f %11.6f ", \
                 num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],                  num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],
                 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++){
Line 1943  void lubksb(double **a, int n, int *indx Line 2003  void lubksb(double **a, int n, int *indx
   }     } 
 }   } 
   
   void pstamp(FILE *fichier)
   {
     fprintf(fichier,"# %s.%s\n#%s\n#%s\n# %s", optionfilefiname,optionfilext,version,fullversion,strstart);
   }
   
 /************ Frequencies ********************/  /************ Frequencies ********************/
 void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[])  void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[])
 {  /* Some frequencies */  {  /* Some frequencies */
Line 1952  void  freqsummary(char fileres[], int ia Line 2017  void  freqsummary(char fileres[], int ia
   double ***freq; /* Frequencies */    double ***freq; /* Frequencies */
   double *pp, **prop;    double *pp, **prop;
   double pos,posprop, k2, dateintsum=0,k2cpt=0;    double pos,posprop, k2, dateintsum=0,k2cpt=0;
   FILE *ficresp;  
   char fileresp[FILENAMELENGTH];    char fileresp[FILENAMELENGTH];
       
   pp=vector(1,nlstate);    pp=vector(1,nlstate);
Line 2017  void  freqsummary(char fileres[], int ia Line 2081  void  freqsummary(char fileres[], int ia
       }        }
                 
       /*      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);*/
 fprintf(ficresp, "#Local time at start: %s", strstart);        pstamp(ficresp);
       if  (cptcovn>0) {        if  (cptcovn>0) {
         fprintf(ficresp, "\n#********** Variable ");           fprintf(ficresp, "\n#********** Variable "); 
         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);          for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
Line 2390  void tricode(int *Tvar, int **nbcode, in Line 2454  void tricode(int *Tvar, int **nbcode, in
   
 /*********** Health Expectancies ****************/  /*********** Health Expectancies ****************/
   
 void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm,double delti[],double **matcov,char strstart[] )  void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] )
   
 {  {
   /* Health expectancies */    /* Health expectancies, no variances */
   int i, j, nhstepm, hstepm, h, nstepm, k, cptj;    int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2;
   double age, agelim, hf;    double age, agelim, hf;
   double ***p3mat,***varhe;    double ***p3mat;
     double eip;
   
     pstamp(ficreseij);
     fprintf(ficreseij,"# (a) Life expectancies by health status at initial age and (b) health expectancies by health status at initial age\n");
     fprintf(ficreseij,"# Age");
     for(i=1; i<=nlstate;i++){
       for(j=1; j<=nlstate;j++){
         fprintf(ficreseij," e%1d%1d ",i,j);
       }
       fprintf(ficreseij," e%1d. ",i);
     }
     fprintf(ficreseij,"\n");
   
     
     if(estepm < stepm){
       printf ("Problem %d lower than %d\n",estepm, stepm);
     }
     else  hstepm=estepm;   
     /* We compute the life expectancy from trapezoids spaced every estepm months
      * This is mainly to measure the difference between two models: for example
      * if stepm=24 months pijx are given only every 2 years and by summing them
      * we are calculating an estimate of the Life Expectancy assuming a linear 
      * progression in between and thus overestimating or underestimating according
      * to the curvature of the survival function. If, for the same date, we 
      * estimate the model with stepm=1 month, we can keep estepm to 24 months
      * to compare the new estimate of Life expectancy with the same linear 
      * hypothesis. A more precise result, taking into account a more precise
      * curvature will be obtained if estepm is as small as stepm. */
   
     /* 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. 
        nhstepm is the number of hstepm from age to agelim 
        nstepm is the number of stepm from age to agelin. 
        Look at hpijx to understand the reason of that which relies in memory size
        and note for a fixed period like estepm months */
     /* 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
        means that if the survival funtion is printed only each two years of age and if
        you sum them up and add 1 year (area under the trapezoids) you won't get the same 
        results. So we changed our mind and took the option of the best precision.
     */
     hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ 
   
     agelim=AGESUP;
     /* nhstepm age range expressed in number of stepm */
     nstepm=(int) rint((agelim-age)*YEARM/stepm); 
     /* Typically if 20 years nstepm = 20*12/6=40 stepm */ 
     /* if (stepm >= YEARM) hstepm=1;*/
     nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   
     for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
       /* Computed by stepm unit matrices, product of hstepm matrices, stored
          in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
       
       hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij);  
       
       hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
       
       printf("%d|",(int)age);fflush(stdout);
       fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
       
       /* Computing expectancies */
       for(i=1; i<=nlstate;i++)
         for(j=1; j<=nlstate;j++)
           for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
             eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf;
             
             /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
   
           }
   
       fprintf(ficreseij,"%3.0f",age );
       for(i=1; i<=nlstate;i++){
         eip=0;
         for(j=1; j<=nlstate;j++){
           eip +=eij[i][j][(int)age];
           fprintf(ficreseij,"%9.4f", eij[i][j][(int)age] );
         }
         fprintf(ficreseij,"%9.4f", eip );
       }
       fprintf(ficreseij,"\n");
       
     }
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
     printf("\n");
     fprintf(ficlog,"\n");
     
   }
   
   void cvevsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] )
   
   {
     /* Covariances of health expectancies eij and of total life expectancies according
      to initial status i, ei. .
     */
     int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji;
     double age, agelim, hf;
     double ***p3matp, ***p3matm, ***varhe;
   double **dnewm,**doldm;    double **dnewm,**doldm;
   double *xp;    double *xp, *xm;
   double **gp, **gm;    double **gp, **gm;
   double ***gradg, ***trgradg;    double ***gradg, ***trgradg;
   int theta;    int theta;
   
     double eip, vip;
   
   varhe=ma3x(1,nlstate*nlstate,1,nlstate*nlstate,(int) bage, (int) fage);    varhe=ma3x(1,nlstate*nlstate,1,nlstate*nlstate,(int) bage, (int) fage);
   xp=vector(1,npar);    xp=vector(1,npar);
     xm=vector(1,npar);
   dnewm=matrix(1,nlstate*nlstate,1,npar);    dnewm=matrix(1,nlstate*nlstate,1,npar);
   doldm=matrix(1,nlstate*nlstate,1,nlstate*nlstate);    doldm=matrix(1,nlstate*nlstate,1,nlstate*nlstate);
       
   fprintf(ficreseij,"# Local time at start: %s", strstart);    pstamp(ficresstdeij);
   fprintf(ficreseij,"# Health expectancies\n");    fprintf(ficresstdeij,"# Health expectancies with standard errors\n");
   fprintf(ficreseij,"# Age");    fprintf(ficresstdeij,"# Age");
   for(i=1; i<=nlstate;i++)    for(i=1; i<=nlstate;i++){
     for(j=1; j<=nlstate;j++)      for(j=1; j<=nlstate;j++)
       fprintf(ficreseij," %1d-%1d (SE)",i,j);        fprintf(ficresstdeij," e%1d%1d (SE)",i,j);
   fprintf(ficreseij,"\n");      fprintf(ficresstdeij," e%1d. ",i);
     }
     fprintf(ficresstdeij,"\n");
   
     pstamp(ficrescveij);
     fprintf(ficrescveij,"# Subdiagonal matrix of covariances of health expectancies by age: cov(eij,ekl)\n");
     fprintf(ficrescveij,"# Age");
     for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++){
         cptj= (j-1)*nlstate+i;
         for(i2=1; i2<=nlstate;i2++)
           for(j2=1; j2<=nlstate;j2++){
             cptj2= (j2-1)*nlstate+i2;
             if(cptj2 <= cptj)
               fprintf(ficrescveij,"  %1d%1d,%1d%1d",i,j,i2,j2);
           }
       }
     fprintf(ficrescveij,"\n");
     
   if(estepm < stepm){    if(estepm < stepm){
     printf ("Problem %d lower than %d\n",estepm, stepm);      printf ("Problem %d lower than %d\n",estepm, stepm);
   }    }
Line 2445  void evsij(char fileres[], double ***eij Line 2628  void evsij(char fileres[], double ***eij
   */    */
   hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */     hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ 
   
     /* If stepm=6 months */
     /* nhstepm age range expressed in number of stepm */
   agelim=AGESUP;    agelim=AGESUP;
   for (age=bage; age<=fage; age ++){ /* If stepm=6 months */    nstepm=(int) rint((agelim-age)*YEARM/stepm); 
     /* nhstepm age range expressed in number of stepm */    /* Typically if 20 years nstepm = 20*12/6=40 stepm */ 
     nstepm=(int) rint((agelim-age)*YEARM/stepm);     /* if (stepm >= YEARM) hstepm=1;*/
     /* Typically if 20 years nstepm = 20*12/6=40 stepm */     nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
     /* if (stepm >= YEARM) hstepm=1;*/    
     nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */    p3matp=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);    p3matm=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
     gradg=ma3x(0,nhstepm,1,npar,1,nlstate*nlstate);    gradg=ma3x(0,nhstepm,1,npar,1,nlstate*nlstate);
     gp=matrix(0,nhstepm,1,nlstate*nlstate);    trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar);
     gm=matrix(0,nhstepm,1,nlstate*nlstate);    gp=matrix(0,nhstepm,1,nlstate*nlstate);
     gm=matrix(0,nhstepm,1,nlstate*nlstate);
   
     for (age=bage; age<=fage; age ++){ 
   
     /* Computed by stepm unit matrices, product of hstepm matrices, stored      /* Computed by stepm unit matrices, product of hstepm matrices, stored
        in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */         in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);    
     
   
     hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */      hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
   
     /* Computing  Variances of health expectancies */      /* Computing  Variances of health expectancies */
       /* Gradient is computed with plus gp and minus gm. Code is duplicated in order to
      for(theta=1; theta <=npar; theta++){         decrease memory allocation */
       for(theta=1; theta <=npar; theta++){
       for(i=1; i<=npar; i++){         for(i=1; i<=npar; i++){ 
         xp[i] = x[i] + (i==theta ?delti[theta]:0);          xp[i] = x[i] + (i==theta ?delti[theta]:0);
           xm[i] = x[i] - (i==theta ?delti[theta]:0);
       }        }
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);          hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij);  
         hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij);  
       
       cptj=0;  
       for(j=1; j<= nlstate; j++){        for(j=1; j<= nlstate; j++){
         for(i=1; i<=nlstate; i++){          for(i=1; i<=nlstate; i++){
           cptj=cptj+1;            for(h=0; h<=nhstepm-1; h++){
           for(h=0, gp[h][cptj]=0.; h<=nhstepm-1; h++){              gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.;
             gp[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;              gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.;
           }            }
         }          }
       }        }
             
              for(ij=1; ij<= nlstate*nlstate; ij++)
       for(i=1; i<=npar; i++)   
         xp[i] = x[i] - (i==theta ?delti[theta]:0);  
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);    
         
       cptj=0;  
       for(j=1; j<= nlstate; j++){  
         for(i=1;i<=nlstate;i++){  
           cptj=cptj+1;  
           for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){  
   
             gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;  
           }  
         }  
       }  
       for(j=1; j<= nlstate*nlstate; j++)  
         for(h=0; h<=nhstepm-1; h++){          for(h=0; h<=nhstepm-1; h++){
           gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];            gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta];
         }          }
      }       }/* End theta */
          
 /* End theta */      
       for(h=0; h<=nhstepm-1; h++)
      trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar);  
   
      for(h=0; h<=nhstepm-1; h++)  
       for(j=1; j<=nlstate*nlstate;j++)        for(j=1; j<=nlstate*nlstate;j++)
         for(theta=1; theta <=npar; theta++)          for(theta=1; theta <=npar; theta++)
           trgradg[h][j][theta]=gradg[h][theta][j];            trgradg[h][j][theta]=gradg[h][theta][j];
            
   
      for(i=1;i<=nlstate*nlstate;i++)       for(ij=1;ij<=nlstate*nlstate;ij++)
       for(j=1;j<=nlstate*nlstate;j++)        for(ji=1;ji<=nlstate*nlstate;ji++)
         varhe[i][j][(int)age] =0.;          varhe[ij][ji][(int)age] =0.;
   
      printf("%d|",(int)age);fflush(stdout);       printf("%d|",(int)age);fflush(stdout);
      fprintf(ficlog,"%d|",(int)age);fflush(ficlog);       fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
Line 2523  void evsij(char fileres[], double ***eij Line 2693  void evsij(char fileres[], double ***eij
       for(k=0;k<=nhstepm-1;k++){        for(k=0;k<=nhstepm-1;k++){
         matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov);          matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov);
         matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]);          matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]);
         for(i=1;i<=nlstate*nlstate;i++)          for(ij=1;ij<=nlstate*nlstate;ij++)
           for(j=1;j<=nlstate*nlstate;j++)            for(ji=1;ji<=nlstate*nlstate;ji++)
             varhe[i][j][(int)age] += doldm[i][j]*hf*hf;              varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;
       }        }
     }      }
     /* Computing expectancies */      /* Computing expectancies */
       hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij);  
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++)        for(j=1; j<=nlstate;j++)
         for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){          for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
           eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf;            eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[i][j][h+1])/2.0*hf;
                       
 /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/            /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
   
         }          }
   
     fprintf(ficreseij,"%3.0f",age );      fprintf(ficresstdeij,"%3.0f",age );
     cptj=0;      for(i=1; i<=nlstate;i++){
         eip=0.;
         vip=0.;
         for(j=1; j<=nlstate;j++){
           eip += eij[i][j][(int)age];
           for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */
             vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age];
           fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) );
         }
         fprintf(ficresstdeij," %9.4f (%.4f)", eip, sqrt(vip));
       }
       fprintf(ficresstdeij,"\n");
   
       fprintf(ficrescveij,"%3.0f",age );
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++){        for(j=1; j<=nlstate;j++){
         cptj++;          cptj= (j-1)*nlstate+i;
         fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) );          for(i2=1; i2<=nlstate;i2++)
             for(j2=1; j2<=nlstate;j2++){
               cptj2= (j2-1)*nlstate+i2;
               if(cptj2 <= cptj)
                 fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]);
             }
       }        }
     fprintf(ficreseij,"\n");      fprintf(ficrescveij,"\n");
         
     free_matrix(gm,0,nhstepm,1,nlstate*nlstate);  
     free_matrix(gp,0,nhstepm,1,nlstate*nlstate);  
     free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*nlstate);  
     free_ma3x(trgradg,0,nhstepm,1,nlstate*nlstate,1,npar);  
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);  
   }    }
     free_matrix(gm,0,nhstepm,1,nlstate*nlstate);
     free_matrix(gp,0,nhstepm,1,nlstate*nlstate);
     free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*nlstate);
     free_ma3x(trgradg,0,nhstepm,1,nlstate*nlstate,1,npar);
     free_ma3x(p3matm,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
     free_ma3x(p3matp,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   printf("\n");    printf("\n");
   fprintf(ficlog,"\n");    fprintf(ficlog,"\n");
   
     free_vector(xm,1,npar);
   free_vector(xp,1,npar);    free_vector(xp,1,npar);
   free_matrix(dnewm,1,nlstate*nlstate,1,npar);    free_matrix(dnewm,1,nlstate*nlstate,1,npar);
   free_matrix(doldm,1,nlstate*nlstate,1,nlstate*nlstate);    free_matrix(doldm,1,nlstate*nlstate,1,nlstate*nlstate);
Line 2616  void varevsij(char optionfilefiname[], d Line 2807  void varevsij(char optionfilefiname[], d
   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);
   fprintf(ficresprobmorprev, "#Local time at start: %s", strstart);    pstamp(ficresprobmorprev);
   fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);    fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
   fprintf(ficresprobmorprev,"# Age cov=%-d",ij);    fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
   for(j=nlstate+1; j<=(nlstate+ndeath);j++){    for(j=nlstate+1; j<=(nlstate+ndeath);j++){
Line 2631  void varevsij(char optionfilefiname[], d Line 2822  void varevsij(char optionfilefiname[], d
   fprintf(fichtm,"\n<br>%s  <br>\n",digitp);    fprintf(fichtm,"\n<br>%s  <br>\n",digitp);
 /*   } */  /*   } */
   varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);    varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
  fprintf(ficresvij, "#Local time at start: %s", strstart);    pstamp(ficresvij);
   fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are the stable prevalence in health states i\n");    fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are ");
     if(popbased==1)
       fprintf(ficresvij,"the age specific prevalence observed in the population i.e cross-sectionally\n in each health state (popbased=1)");
     else
       fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n");
   fprintf(ficresvij,"# Age");    fprintf(ficresvij,"# Age");
   for(i=1; i<=nlstate;i++)    for(i=1; i<=nlstate;i++)
     for(j=1; j<=nlstate;j++)      for(j=1; j<=nlstate;j++)
       fprintf(ficresvij," Cov(e%1d, e%1d)",i,j);        fprintf(ficresvij," Cov(e.%1d, e.%1d)",i,j);
   fprintf(ficresvij,"\n");    fprintf(ficresvij,"\n");
   
   xp=vector(1,npar);    xp=vector(1,npar);
Line 2877  void varprevlim(char fileres[], double * Line 3072  void varprevlim(char fileres[], double *
   double **gradg, **trgradg;    double **gradg, **trgradg;
   double age,agelim;    double age,agelim;
   int theta;    int theta;
   fprintf(ficresvpl, "#Local time at start: %s", strstart);     
   fprintf(ficresvpl,"# Standard deviation of stable prevalences \n");    pstamp(ficresvpl);
     fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n");
   fprintf(ficresvpl,"# Age");    fprintf(ficresvpl,"# Age");
   for(i=1; i<=nlstate;i++)    for(i=1; i<=nlstate;i++)
       fprintf(ficresvpl," %1d-%1d",i,i);        fprintf(ficresvpl," %1d-%1d",i,i);
Line 2992  void varprob(char optionfilefiname[], do Line 3188  void varprob(char optionfilefiname[], do
   fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);    fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
   printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);    printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
   fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);    fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
   fprintf(ficresprob, "#Local time at start: %s", strstart);    pstamp(ficresprob);
   fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");    fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
   fprintf(ficresprob,"# Age");    fprintf(ficresprob,"# Age");
   fprintf(ficresprobcov, "#Local time at start: %s", strstart);    pstamp(ficresprobcov);
   fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");    fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");
   fprintf(ficresprobcov,"# Age");    fprintf(ficresprobcov,"# Age");
   fprintf(ficresprobcor, "#Local time at start: %s", strstart);    pstamp(ficresprobcor);
   fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");    fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");
   fprintf(ficresprobcov,"# Age");    fprintf(ficresprobcor,"# Age");
   
   
   for(i=1; i<=nlstate;i++)    for(i=1; i<=nlstate;i++)
Line 3281  void printinghtml(char fileres[], char t Line 3477  void printinghtml(char fileres[], char t
  - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",   - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
            stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij"));             stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij"));
    fprintf(fichtm,"\     fprintf(fichtm,"\
  - 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(fileres,"pl"),subdirf2(fileres,"pl"));             subdirf2(fileres,"pl"),subdirf2(fileres,"pl"));
    fprintf(fichtm,"\     fprintf(fichtm,"\
  - Life expectancies by age and initial health status (estepm=%2d months) WRONG LINK (to be made): \   - (a) Life expectancies by health status at initial age, (b) health expectancies by health status at initial age:  ei., eij . If one or more covariate 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</li>",  
            estepm,subdirf2(fileres,"le"),subdirf2(fileres,"le"));  
    fprintf(fichtm,"\  
  - Health expectancies by age and initial health status with standard errors (estepm=%2d months): \  
    <a href=\"%s\">%s</a> <br>\n</li>",     <a href=\"%s\">%s</a> <br>\n</li>",
            estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));             estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));
    fprintf(fichtm,"\  
  - Variances and covariances of health expectancies by age and initial health status (estepm=%2d months) TO BE MADE: \  
    <a href=\"%s\">%s</a> <br>\n</li>",  
            estepm,subdirf2(fileres,"vch"),subdirf2(fileres,"vch"));  
   
   
 fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");  fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
Line 3313  fprintf(fichtm," \n<ul><li><b>Graphs</b> Line 3501  fprintf(fichtm," \n<ul><li><b>Graphs</b>
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }       }
      /* Pij */       /* Pij */
      fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s%d1.png<br> \       fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s%d1.png\">%s%d1.png</a><br> \
 <img src=\"%s%d1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);       <img src=\"%s%d1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);     
      /* Quasi-incidences */       /* Quasi-incidences */
      fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\       fprintf(fichtm,"<br>- Pij 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: %s%d2.png<br> \   before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: <a href=\"%s%d2.png\">%s%d2.png</a><br> \
 <img src=\"%s%d2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);   <img src=\"%s%d2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); 
        /* 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>- Stable prevalence in each health state : p%s%d%d.png<br> \           fprintf(fichtm,"<br>- Period (stable) prevalence in each health state : <a href=\"%s%d%d.png\">%s%d%d.png</a><br> \
 <img src=\"%s%d%d.png\">",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);  <img src=\"%s%d%d.png\">",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
        }         }
      for(cpt=1; cpt<=nlstate;cpt++) {       for(cpt=1; cpt<=nlstate;cpt++) {
         fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): %s%d%d.png <br> \          fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \
 <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);  <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
      }       }
    } /* end i1 */     } /* end i1 */
  }/* End k1 */   }/* End k1 */
Line 3347  fprintf(fichtm," \n<ul><li><b>Graphs</b> Line 3535  fprintf(fichtm," \n<ul><li><b>Graphs</b>
  - Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",   - Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
          subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor"));           subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor"));
  fprintf(fichtm,"\   fprintf(fichtm,"\
  - Variances and covariances of health expectancies by age (estepm=%d months): <a href=\"%s\">%s</a><br>\n",   - Variances and covariances of health expectancies by age and <b>initial health status</b> (cov(e<sup>ij</sup>,e<sup>kl</sup>)(estepm=%2d months): \
      <a href=\"%s\">%s</a> <br>\n</li>",
              estepm,subdirf2(fileres,"cve"),subdirf2(fileres,"cve"));
    fprintf(fichtm,"\
    - (a) Health expectancies by health status at initial age (e<sup>ij</sup>) and standard errors (in parentheses) (b) life expectancies and standard errors (e<sup>i.</sup>=e<sup>i1</sup>+e<sup>i2</sup>+...)(estepm=%2d months): \
      <a href=\"%s\">%s</a> <br>\n</li>",
              estepm,subdirf2(fileres,"stde"),subdirf2(fileres,"stde"));
    fprintf(fichtm,"\
    - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences (i.e population based) (estepm=%d months): <a href=\"%s\">%s</a><br>\n",
          estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"));           estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"));
  fprintf(fichtm,"\   fprintf(fichtm,"\
  - Life and health expectancies with their standard errors: <a href=\"%s\">%s</a> <br>\n",   - Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors: <a href=\"%s\">%s</a> <br>\n",
          subdirf2(fileres,"t"),subdirf2(fileres,"t"));           subdirf2(fileres,"t"),subdirf2(fileres,"t"));
  fprintf(fichtm,"\   fprintf(fichtm,"\
  - Standard deviation of stable prevalences: <a href=\"%s\">%s</a> <br>\n",\   - Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\
          subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl"));           subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl"));
   
 /*  if(popforecast==1) fprintf(fichtm,"\n */  /*  if(popforecast==1) fprintf(fichtm,"\n */
Line 3425  plot [%.f:%.f] \"%s\" every :::%d::%d u Line 3621  plot [%.f:%.f] \"%s\" every :::%d::%d u
        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\"Stable prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);       fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
      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 3473  plot [%.f:%.f] \"%s\" every :::%d::%d u Line 3669  plot [%.f:%.f] \"%s\" every :::%d::%d u
       
   for (k1=1; k1<= m ; k1 ++) {     for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<= nlstate ; cpt ++) {      for (cpt=1; cpt<= nlstate ; cpt ++) {
       k=2+nlstate*(2*cpt-2);        /*       k=2+nlstate*(2*cpt-2); */
         k=2+(nlstate+1)*(cpt-1);
       fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1);        fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1);
       fprintf(ficgp,"set ter png small\n\        fprintf(ficgp,"set ter png small\n\
 set size 0.65,0.65\n\  set size 0.65,0.65\n\
Line 3487  plot [%.f:%.f] \"%s\" every :::%d::%d u Line 3684  plot [%.f:%.f] \"%s\" every :::%d::%d u
                   
       */        */
       for (i=1; i< nlstate ; i ++) {        for (i=1; i< nlstate ; i ++) {
         fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);          fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+i,cpt,i+1);
           /*      fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
                   
       }         } 
         fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+nlstate,cpt);
     }      }
   }    }
       
Line 4130  void printinggnuplotmort(char fileres[], Line 4329  void printinggnuplotmort(char fileres[],
   
   
   
   
 /***********************************************/  /***********************************************/
 /**************** Main Program *****************/  /**************** Main Program *****************/
 /***********************************************/  /***********************************************/
Line 4198  int main(int argc, char *argv[]) Line 4398  int main(int argc, char *argv[])
   char z[1]="c", occ;    char z[1]="c", occ;
   
   char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];    char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];
   char strstart[80], *strt, strtend[80];    char  *strt, strtend[80];
   char *stratrunc;    char *stratrunc;
   int lstra;    int lstra;
   
Line 4249  int main(int argc, char *argv[]) Line 4449  int main(int argc, char *argv[])
       while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0');        while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0');
       printf("val= |%s| pathr=%s\n",val,pathr);        printf("val= |%s| pathr=%s\n",val,pathr);
       strcpy (pathtot, val);        strcpy (pathtot, val);
       if(pathr[0] == '\0') break; /* Un peu sale */        if(pathr[0] == '\0') break; /* Dirty */
     }      }
   }    }
   else{    else{
Line 4268  int main(int argc, char *argv[]) Line 4468  int main(int argc, char *argv[])
   /* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */    /* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */
   split(pathtot,path,optionfile,optionfilext,optionfilefiname);    split(pathtot,path,optionfile,optionfilext,optionfilefiname);
   printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);    printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
   chdir(path);    chdir(path); /* Can be a relative path */
     if(getcwd(pathcd,MAXLINE) > 0) /* So pathcd is the full path */
       printf("Current directory %s!\n",pathcd);
   strcpy(command,"mkdir ");    strcpy(command,"mkdir ");
   strcat(command,optionfilefiname);    strcat(command,optionfilefiname);
   if((outcmd=system(command)) != 0){    if((outcmd=system(command)) != 0){
Line 4377  int main(int argc, char *argv[]) Line 4579  int main(int argc, char *argv[])
     free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);       free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); 
     fclose (ficparo);      fclose (ficparo);
     fclose (ficlog);      fclose (ficlog);
       goto end;
     exit(0);      exit(0);
   }    }
   else if(mle==-3) {    else if(mle==-3) {
Line 4407  int main(int argc, char *argv[]) Line 4610  int main(int argc, char *argv[])
         j++;          j++;
         fscanf(ficpar,"%1d%1d",&i1,&j1);          fscanf(ficpar,"%1d%1d",&i1,&j1);
         if ((i1 != i) && (j1 != j)){          if ((i1 != i) && (j1 != j)){
           printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);            printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \
   It might be a problem of design; if ncovcol and the model are correct\n \
   run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1);
           exit(1);            exit(1);
         }          }
         fprintf(ficparo,"%1d%1d",i1,j1);          fprintf(ficparo,"%1d%1d",i1,j1);
Line 4534  int main(int argc, char *argv[]) Line 4739  int main(int argc, char *argv[])
   
   /*-------- data file ----------*/    /*-------- data file ----------*/
   if((fic=fopen(datafile,"r"))==NULL)    {    if((fic=fopen(datafile,"r"))==NULL)    {
     printf("Problem with datafile: %s\n", datafile);goto end;      printf("Problem while opening datafile: %s\n", datafile);goto end;
     fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end;      fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);goto end;
   }    }
   
   n= lastobs;    n= lastobs;
Line 4631  int main(int argc, char *argv[]) Line 4836  int main(int argc, char *argv[])
           
     cutv(stra, strb,line,' ');       cutv(stra, strb,line,' '); 
     errno=0;      errno=0;
     lval=strtol(strb,&endptr,10);       dval=strtod(strb,&endptr); 
     if( strb[0]=='\0' || (*endptr != '\0')){      if( strb[0]=='\0' || (*endptr != '\0')){
       printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a weight.  Exiting.\n",lval, i,line,linei);        printf("Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight.  Exiting.\n",dval, i,line,linei);
       exit(1);        exit(1);
     }      }
     weight[i]=(double)(lval);       weight[i]=dval; 
     strcpy(line,stra);      strcpy(line,stra);
           
     for (j=ncovcol;j>=1;j--){      for (j=ncovcol;j>=1;j--){
Line 4648  int main(int argc, char *argv[]) Line 4853  int main(int argc, char *argv[])
         exit(1);          exit(1);
       }        }
       if(lval <-1 || lval >1){        if(lval <-1 || lval >1){
         printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a value of the %d covar (meaning 0 for the reference or 1. IMaCh does not build design variables, do it your self).  Exiting.\n",lval,linei, i,line,j);          printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\n \
    Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
    for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
    For example, for multinomial values like 1, 2 and 3,\n \
    build V1=0 V2=0 for the reference value (1),\n \
           V1=1 V2=0 for (2) \n \
    and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
    output of IMaCh is often meaningless.\n \
    Exiting.\n",lval,linei, i,line,j);
         exit(1);          exit(1);
       }        }
       covar[j][i]=(double)(lval);        covar[j][i]=(double)(lval);
Line 4944  int main(int argc, char *argv[]) Line 5157  int main(int argc, char *argv[])
     printf("Problem with %s \n",optionfilehtmcov), exit(0);      printf("Problem with %s \n",optionfilehtmcov), exit(0);
   }    }
   else{    else{
   fprintf(fichtmcov,"<body>\n<title>IMaCh Cov %s</title>\n <font size=\"2\">%s <br> %s</font> \    fprintf(fichtmcov,"<html><head>\n<title>IMaCh Cov %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n",\  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n",\
           fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);            optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }    }
   
   fprintf(fichtm,"<body>\n<title>IMaCh %s</title>\n <font size=\"2\">%s <br> %s</font> \    fprintf(fichtm,"<html><head>\n<title>IMaCh %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n\  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n\
 \n\  \n\
 <hr  size=\"2\" color=\"#EC5E5E\">\  <hr  size=\"2\" color=\"#EC5E5E\">\
  <ul><li><h4>Parameter files</h4>\n\   <ul><li><h4>Parameter files</h4>\n\
    - Parameter file: <a href=\"%s.%s\">%s.%s</a><br>\n\
  - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n\   - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n\
  - Log file of the run: <a href=\"%s\">%s</a><br>\n\   - Log file of the run: <a href=\"%s\">%s</a><br>\n\
  - Gnuplot file name: <a href=\"%s\">%s</a><br>\n\   - Gnuplot file name: <a href=\"%s\">%s</a><br>\n\
  - Date and time at start: %s</ul>\n",\   - Date and time at start: %s</ul>\n",\
           fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\            optionfilehtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\
             optionfilefiname,optionfilext,optionfilefiname,optionfilext,\
           fileres,fileres,\            fileres,fileres,\
           filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);            filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);
   fflush(fichtm);    fflush(fichtm);
Line 5099  Interval (in months) between two waves: Line 5314  Interval (in months) between two waves:
       printf("%d %.0lf %lf %.0lf %.0lf %.0lf %lf\n",k,lsurv[k],p[1]*exp(p[2]*(k-agegomp)),(p[1]*exp(p[2]*(k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]);        printf("%d %.0lf %lf %.0lf %.0lf %.0lf %lf\n",k,lsurv[k],p[1]*exp(p[2]*(k-agegomp)),(p[1]*exp(p[2]*(k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]);
           
           
     replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */      replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */
     printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);      printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
           
     printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \      printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \
Line 5336  Interval (in months) between two waves: Line 5551  Interval (in months) between two waves:
     /*  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);*/      /*  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);*/
     /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/      /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
           
     replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */      replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */
     printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);      printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
           
     printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\      printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\
Line 5358  Interval (in months) between two waves: Line 5573  Interval (in months) between two waves:
     fclose(ficres);      fclose(ficres);
   
   
     /*--------------- Prevalence limit  (stable prevalence) --------------*/      /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
       
     strcpy(filerespl,"pl");      strcpy(filerespl,"pl");
     strcat(filerespl,fileres);      strcat(filerespl,fileres);
     if((ficrespl=fopen(filerespl,"w"))==NULL) {      if((ficrespl=fopen(filerespl,"w"))==NULL) {
       printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end;        printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end;
       fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end;        fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end;
     }      }
     printf("Computing stable prevalence: result on file '%s' \n", filerespl);      printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl);
     fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl);      fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl);
     fprintf(ficrespl, "#Local time at start: %s", strstart);      pstamp(ficrespl);
     fprintf(ficrespl,"#Stable prevalence \n");      fprintf(ficrespl,"# Period (stable) prevalence \n");
     fprintf(ficrespl,"#Age ");      fprintf(ficrespl,"#Age ");
     for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);      for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
     fprintf(ficrespl,"\n");      fprintf(ficrespl,"\n");
Line 5429  Interval (in months) between two waves: Line 5644  Interval (in months) between two waves:
     hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */       hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ 
   
     /* hstepm=1;   aff par mois*/      /* hstepm=1;   aff par mois*/
     fprintf(ficrespij, "#Local time at start: %s", strstart);      pstamp(ficrespij);
     fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");      fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
     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 5499  Interval (in months) between two waves: Line 5714  Interval (in months) between two waves:
       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 LEs with variances: file '%s' \n", filerest);       printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); 
     fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest);       fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); 
   
   
     strcpy(filerese,"e");      strcpy(filerese,"e");
Line 5512  Interval (in months) between two waves: Line 5727  Interval (in months) between two waves:
     printf("Computing Health Expectancies: result on file '%s' \n", filerese);      printf("Computing Health Expectancies: result on file '%s' \n", filerese);
     fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);      fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
   
       strcpy(fileresstde,"stde");
       strcat(fileresstde,fileres);
       if((ficresstdeij=fopen(fileresstde,"w"))==NULL) {
         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);
       }
       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);
   
       strcpy(filerescve,"cve");
       strcat(filerescve,fileres);
       if((ficrescveij=fopen(filerescve,"w"))==NULL) {
         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);
       }
       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);
   
     strcpy(fileresv,"v");      strcpy(fileresv,"v");
     strcat(fileresv,fileres);      strcat(fileresv,fileres);
     if((ficresvij=fopen(fileresv,"w"))==NULL) {      if((ficresvij=fopen(fileresv,"w"))==NULL) {
Line 5544  Interval (in months) between two waves: Line 5777  Interval (in months) between two waves:
         fprintf(ficrest,"******\n");          fprintf(ficrest,"******\n");
   
         fprintf(ficreseij,"\n#****** ");          fprintf(ficreseij,"\n#****** ");
         for(j=1;j<=cptcoveff;j++)           fprintf(ficresstdeij,"\n#****** ");
           fprintf(ficrescveij,"\n#****** ");
           for(j=1;j<=cptcoveff;j++) {
           fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);            fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
             fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
             fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
           }
         fprintf(ficreseij,"******\n");          fprintf(ficreseij,"******\n");
           fprintf(ficresstdeij,"******\n");
           fprintf(ficrescveij,"******\n");
   
         fprintf(ficresvij,"\n#****** ");          fprintf(ficresvij,"\n#****** ");
         for(j=1;j<=cptcoveff;j++)           for(j=1;j<=cptcoveff;j++) 
Line 5555  Interval (in months) between two waves: Line 5795  Interval (in months) between two waves:
   
         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;
         evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);            evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);  
           cvevsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);  
     
         vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);          vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
         oldm=oldms;savm=savms;          oldm=oldms;savm=savms;
Line 5564  Interval (in months) between two waves: Line 5805  Interval (in months) between two waves:
           varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav, strstart);            varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav, strstart);
         }          }
   
         fprintf(ficrest, "#Local time at start: %s", strstart);          pstamp(ficrest);
         fprintf(ficrest,"#Total LEs with variances: e.. (std) ");          fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# Age ( e.. (std) ");
         for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);          for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
         fprintf(ficrest,"\n");          fprintf(ficrest,"\n");
   
Line 5613  Interval (in months) between two waves: Line 5854  Interval (in months) between two waves:
     free_ivector(cod,1,n);      free_ivector(cod,1,n);
     free_ivector(tab,1,NCOVMAX);      free_ivector(tab,1,NCOVMAX);
     fclose(ficreseij);      fclose(ficreseij);
       fclose(ficresstdeij);
       fclose(ficrescveij);
     fclose(ficresvij);      fclose(ficresvij);
     fclose(ficrest);      fclose(ficrest);
     fclose(ficpar);      fclose(ficpar);
       
     /*------- Variance of stable prevalence------*/         /*------- Variance of period (stable) prevalence------*/   
   
     strcpy(fileresvpl,"vpl");      strcpy(fileresvpl,"vpl");
     strcat(fileresvpl,fileres);      strcat(fileresvpl,fileres);
     if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {      if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
       printf("Problem with variance of 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 stable prevalence: file '%s' \n", fileresvpl);      printf("Computing Variance-covariance of period (stable) prevalence: file '%s' \n", fileresvpl);
   
     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 5696  Interval (in months) between two waves: Line 5939  Interval (in months) between two waves:
   fprintf(ficlog,"Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec);    fprintf(ficlog,"Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec);
   /*  printf("Total time was %d uSec.\n", total_usecs);*/    /*  printf("Total time was %d uSec.\n", total_usecs);*/
 /*   if(fileappend(fichtm,optionfilehtm)){ */  /*   if(fileappend(fichtm,optionfilehtm)){ */
   fprintf(fichtm,"<br>Local time at start %s<br>Local time at end   %s<br>",strstart, strtend);    fprintf(fichtm,"<br>Local time at start %s<br>Local time at end   %s<br>\n</body></html>",strstart, strtend);
   fclose(fichtm);    fclose(fichtm);
     fprintf(fichtmcov,"<br>Local time at start %s<br>Local time at end   %s<br>\n</body></html>",strstart, strtend);
   fclose(fichtmcov);    fclose(fichtmcov);
   fclose(ficgp);    fclose(ficgp);
   fclose(ficlog);    fclose(ficlog);
   /*------ End -----------*/    /*------ End -----------*/
   
   chdir(path);  
      printf("Before Current directory %s!\n",pathcd);
      if(chdir(pathcd) != 0)
       printf("Can't move to directory %s!\n",path);
     if(getcwd(pathcd,MAXLINE) > 0)
       printf("Current directory %s!\n",pathcd);
   /*strcat(plotcmd,CHARSEPARATOR);*/    /*strcat(plotcmd,CHARSEPARATOR);*/
   sprintf(plotcmd,"gnuplot");    sprintf(plotcmd,"gnuplot");
 #ifndef UNIX  #ifndef UNIX
Removed from v.1.116  
changed lines
  Added in v.1.123

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