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

Diff for /imach/src/imach.c between versions 1.41.2.2 and 1.48

version 1.41.2.2, 2003/06/13 07:45:28 version 1.48, 2002/06/10 13:12:49
Line 56 Line 56
 #include <unistd.h>  #include <unistd.h>
   
 #define MAXLINE 256  #define MAXLINE 256
 #define GNUPLOTPROGRAM "wgnuplot"  #define GNUPLOTPROGRAM "gnuplot"
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/  /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
 #define FILENAMELENGTH 80  #define FILENAMELENGTH 80
 /*#define DEBUG*/  /*#define DEBUG*/
   #define windows
 /*#define windows*/  
 #define GLOCK_ERROR_NOPATH              -1      /* empty path */  #define GLOCK_ERROR_NOPATH              -1      /* empty path */
 #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */  #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */
   
Line 76 Line 75
 #define YEARM 12. /* Number of months per year */  #define YEARM 12. /* Number of months per year */
 #define AGESUP 130  #define AGESUP 130
 #define AGEBASE 40  #define AGEBASE 40
   #ifdef windows
   #define DIRSEPARATOR '\\'
   #else
   #define DIRSEPARATOR '/'
   #endif
   
   char version[80]="Imach version 0.8h, May 2002, INED-EUROREVES ";
 int erreur; /* Error number */  int erreur; /* Error number */
 int nvar;  int nvar;
 int cptcovn, cptcovage=0, cptcoveff=0,cptcov;  int cptcovn, cptcovage=0, cptcoveff=0,cptcov;
Line 97  double jmean; /* Mean space between 2 wa Line 101  double jmean; /* Mean space between 2 wa
 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,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
 FILE *ficgp,*ficresprob,*ficpop;  FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;
   FILE *fichtm; /* Html File */
 FILE *ficreseij;  FILE *ficreseij;
   char filerese[FILENAMELENGTH];  char filerese[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 optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];
   char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
   
   char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
   
   char filerest[FILENAMELENGTH];
   char fileregp[FILENAMELENGTH];
   char popfile[FILENAMELENGTH];
   
   char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];
   
 #define NR_END 1  #define NR_END 1
 #define FREE_ARG char*  #define FREE_ARG char*
Line 162  static int split( char *path, char *dirc Line 178  static int split( char *path, char *dirc
   
    l1 = strlen( path );                 /* length of path */     l1 = strlen( path );                 /* length of path */
    if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );     if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
 #ifdef windows     s = strrchr( path,  DIRSEPARATOR );          /* find last / */
    s = strrchr( path, '\\' );           /* find last / */  
 #else  
    s = strrchr( path, '/' );            /* find last / */  
 #endif  
    if ( s == NULL ) {                   /* no directory, so use current */     if ( s == NULL ) {                   /* no directory, so use current */
 #if     defined(__bsd__)                /* get current working directory */  #if     defined(__bsd__)                /* get current working directory */
       extern char       *getwd( );        extern char       *getwd( );
Line 870  double func( double *x) Line 882  double func( double *x)
   double **out;    double **out;
   double sw; /* Sum of weights */    double sw; /* Sum of weights */
   double lli; /* Individual log likelihood */    double lli; /* Individual log likelihood */
   int s1, s2;  
   long ipmx;    long ipmx;
   /*extern weight */    /*extern weight */
   /* We are differentiating ll according to initial status */    /* We are differentiating ll according to initial status */
Line 885  double func( double *x) Line 896  double func( double *x)
     for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];      for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
     for(mi=1; mi<= wav[i]-1; mi++){      for(mi=1; mi<= wav[i]-1; mi++){
       for (ii=1;ii<=nlstate+ndeath;ii++)        for (ii=1;ii<=nlstate+ndeath;ii++)
         for (j=1;j<=nlstate+ndeath;j++){          for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0);
           oldm[ii][j]=(ii==j ? 1.0 : 0.0);  
           savm[ii][j]=(ii==j ? 1.0 : 0.0);  
         }  
       for(d=0; d<dh[mi][i]; d++){        for(d=0; d<dh[mi][i]; d++){
         newm=savm;          newm=savm;
         cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;          cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
Line 904  double func( double *x) Line 912  double func( double *x)
                 
       } /* end mult */        } /* end mult */
             
       s1=s[mw[mi][i]][i];        lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);
       s2=s[mw[mi+1][i]][i];        /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/
       if( s2 > nlstate){  
         /* i.e. if s2 is a death state and if the date of death is known then the contribution  
            to the likelihood is the probability to die between last step unit time and current  
            step unit time, which is also the differences between probability to die before dh  
            and probability to die before dh-stepm .  
            In version up to 0.92 likelihood was computed  
            as if date of death was unknown. Death was treated as any other  
            health state: the date of the interview describes the actual state  
            and not the date of a change in health state. The former idea was  
            to consider that at each interview the state was recorded  
            (healthy, disable or death) and IMaCh was corrected; but when we  
            introduced the exact date of death then we should have modified  
            the contribution of an exact death to the likelihood. This new  
            contribution is smaller and very dependent of the step unit  
            stepm. It is no more the probability to die between last interview  
            and month of death but the probability to survive from last  
            interview up to one month before death multiplied by the  
            probability to die within a month. Thanks to Chris  
            Jackson for correcting this bug.  Former versions increased  
            mortality artificially. The bad side is that we add another loop  
            which slows down the processing. The difference can be up to 10%  
            lower mortality.  
         */  
         lli=log(out[s1][s2] - savm[s1][s2]);  
       }else{  
         lli=log(out[s1][s2]); /* or     lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); */  
         /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/  
       }  
       ipmx +=1;        ipmx +=1;
       sw += weight[i];        sw += weight[i];
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;        ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
       /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d lli=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],lli,weight[i],out[s1][s2],savm[s1][s2]);*/  
     } /* end of wave */      } /* end of wave */
   } /* end of individual */    } /* end of individual */
   
   for(k=1,l=0.; k<=nlstate; k++) l += ll[k];    for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
   /* printf("l1=%f l2=%f ",ll[1],ll[2]); */    /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
   l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */    l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
   /*exit(0);*/  
   return -l;    return -l;
 }  }
   
Line 1365  void prevalence(int agemin, float agemax Line 1343  void prevalence(int agemin, float agemax
   j=cptcoveff;    j=cptcoveff;
   if (cptcovn<1) {j=1;ncodemax[1]=1;}    if (cptcovn<1) {j=1;ncodemax[1]=1;}
     
  for(k1=1; k1<=j;k1++){    for(k1=1; k1<=j;k1++){
     for(i1=1; i1<=ncodemax[k1];i1++){      for(i1=1; i1<=ncodemax[k1];i1++){
       j1++;        j1++;
         
       for (i=-1; i<=nlstate+ndeath; i++)          for (i=-1; i<=nlstate+ndeath; i++)  
         for (jk=-1; jk<=nlstate+ndeath; jk++)            for (jk=-1; jk<=nlstate+ndeath; jk++)  
           for(m=agemin; m <= agemax+3; m++)            for(m=agemin; m <= agemax+3; m++)
Line 1387  void prevalence(int agemin, float agemax Line 1365  void prevalence(int agemin, float agemax
             if ((k2>=dateprev1) && (k2<=dateprev2)) {              if ((k2>=dateprev1) && (k2<=dateprev2)) {
               if(agev[m][i]==0) agev[m][i]=agemax+1;                if(agev[m][i]==0) agev[m][i]=agemax+1;
               if(agev[m][i]==1) agev[m][i]=agemax+2;                if(agev[m][i]==1) agev[m][i]=agemax+2;
               if (m<lastpass)                if (m<lastpass) {
                 if (calagedate>0) freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];                  if (calagedate>0)
               else                    freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
                freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];                  else
                freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i];                    freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
                   freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i];
                 }
             }              }
           }            }
         }          }
       }        }
         for(i=(int)agemin; i <= (int)agemax+3; i++){        for(i=(int)agemin; i <= (int)agemax+3; i++){
           for(jk=1; jk <=nlstate ; jk++){          for(jk=1; jk <=nlstate ; jk++){
             for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)            for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
               pp[jk] += freq[jk][m][i];              pp[jk] += freq[jk][m][i];
           }          }
           for(jk=1; jk <=nlstate ; jk++){          for(jk=1; jk <=nlstate ; jk++){
             for(m=-1, pos=0; m <=0 ; m++)            for(m=-1, pos=0; m <=0 ; m++)
             pos += freq[jk][m][i];              pos += freq[jk][m][i];
         }          }
                 
          for(jk=1; jk <=nlstate ; jk++){          for(jk=1; jk <=nlstate ; jk++){
            for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)            for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
              pp[jk] += freq[jk][m][i];              pp[jk] += freq[jk][m][i];
          }          }
                   
          for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];          for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];
          
          for(jk=1; jk <=nlstate ; jk++){                    for(jk=1; jk <=nlstate ; jk++){    
            if( i <= (int) agemax){            if( i <= (int) agemax){
              if(pos>=1.e-5){              if(pos>=1.e-5){
                probs[i][jk][j1]= pp[jk]/pos;                probs[i][jk][j1]= pp[jk]/pos;
              }              }
            }            }
          }  
            
         }          }
          
         }
     }      }
   }    }
   
Line 1674  void evsij(char fileres[], double ***eij Line 1654  void evsij(char fileres[], double ***eij
           }            }
         }          }
       }        }
        
      
   
       for(j=1; j<= nlstate*2; j++)        for(j=1; j<= nlstate*2; 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][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
         }          }
   
      }       }
         
 /* End theta */  /* End theta */
Line 1691  void evsij(char fileres[], double ***eij Line 1667  void evsij(char fileres[], double ***eij
      for(h=0; h<=nhstepm-1; h++)       for(h=0; h<=nhstepm-1; h++)
       for(j=1; j<=nlstate*2;j++)        for(j=1; j<=nlstate*2;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*2;i++)       for(i=1;i<=nlstate*2;i++)
       for(j=1;j<=nlstate*2;j++)        for(j=1;j<=nlstate*2;j++)
         varhe[i][j][(int)age] =0.;          varhe[i][j][(int)age] =0.;
   
     for(h=0;h<=nhstepm-1;h++){       printf("%d|",(int)age);fflush(stdout);
        for(h=0;h<=nhstepm-1;h++){
       for(k=0;k<=nhstepm-1;k++){        for(k=0;k<=nhstepm-1;k++){
         matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);          matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);
         matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]);          matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]);
Line 1707  void evsij(char fileres[], double ***eij Line 1684  void evsij(char fileres[], double ***eij
             varhe[i][j][(int)age] += doldm[i][j]*hf*hf;              varhe[i][j][(int)age] += doldm[i][j]*hf*hf;
       }        }
     }      }
   
        
     /* Computing expectancies */      /* Computing expectancies */
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++)        for(j=1; j<=nlstate;j++)
Line 1734  void evsij(char fileres[], double ***eij Line 1709  void evsij(char fileres[], double ***eij
     free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar);      free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar);
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   }    }
     printf("\n");
   
   free_vector(xp,1,npar);    free_vector(xp,1,npar);
   free_matrix(dnewm,1,nlstate*2,1,npar);    free_matrix(dnewm,1,nlstate*2,1,npar);
   free_matrix(doldm,1,nlstate*2,1,nlstate*2);    free_matrix(doldm,1,nlstate*2,1,nlstate*2);
Line 1756  void varevsij(char fileres[], double *** Line 1733  void varevsij(char fileres[], double ***
   double age,agelim, hf;    double age,agelim, hf;
   int theta;    int theta;
   
    fprintf(ficresvij,"# Covariances of life expectancies\n");    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,"# 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++)
Line 1891  void varprevlim(char fileres[], double * Line 1868  void varprevlim(char fileres[], double *
   double age,agelim;    double age,agelim;
   int theta;    int theta;
         
   fprintf(ficresvpl,"# Standard deviation of prevalences limit\n");    fprintf(ficresvpl,"# Standard deviation of prevalence's limit\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 1960  void varprevlim(char fileres[], double * Line 1937  void varprevlim(char fileres[], double *
 }  }
   
 /************ Variance of one-step probabilities  ******************/  /************ Variance of one-step probabilities  ******************/
 void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)  void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)
 {  {
   int i, j, i1, k1, j1, z1;    int i, j,  i1, k1, l1;
   int k=0, cptcode;    int k2, l2, j1,  z1;
     int k=0,l, cptcode;
     int first=1;
     double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2;
   double **dnewm,**doldm;    double **dnewm,**doldm;
   double *xp;    double *xp;
   double *gp, *gm;    double *gp, *gm;
   double **gradg, **trgradg;    double **gradg, **trgradg;
     double **mu;
   double age,agelim, cov[NCOVMAX];    double age,agelim, cov[NCOVMAX];
     double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */
   int theta;    int theta;
   char fileresprob[FILENAMELENGTH];    char fileresprob[FILENAMELENGTH];
     char fileresprobcov[FILENAMELENGTH];
     char fileresprobcor[FILENAMELENGTH];
   
     double ***varpij;
   
   strcpy(fileresprob,"prob");    strcpy(fileresprob,"prob");
   strcat(fileresprob,fileres);    strcat(fileresprob,fileres);
   if((ficresprob=fopen(fileresprob,"w"))==NULL) {    if((ficresprob=fopen(fileresprob,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprob);      printf("Problem with resultfile: %s\n", fileresprob);
   }    }
     strcpy(fileresprobcov,"probcov");
     strcat(fileresprobcov,fileres);
     if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", fileresprobcov);
     }
     strcpy(fileresprobcor,"probcor");
     strcat(fileresprobcor,fileres);
     if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", fileresprobcor);
     }
   printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);    printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
     printf("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);
     
 fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n");    fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
   fprintf(ficresprob,"# Age");    fprintf(ficresprob,"# Age");
   for(i=1; i<=nlstate;i++)    fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");
     for(j=1; j<=(nlstate+ndeath);j++)    fprintf(ficresprobcov,"# Age");
       fprintf(ficresprob," p%1d-%1d (SE)",i,j);    fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");
     fprintf(ficresprobcov,"# Age");
   
   
     for(i=1; i<=nlstate;i++)
       for(j=1; j<=(nlstate+ndeath);j++){
         fprintf(ficresprob," p%1d-%1d (SE)",i,j);
         fprintf(ficresprobcov," p%1d-%1d ",i,j);
         fprintf(ficresprobcor," p%1d-%1d ",i,j);
       }  
   fprintf(ficresprob,"\n");    fprintf(ficresprob,"\n");
     fprintf(ficresprobcov,"\n");
     fprintf(ficresprobcor,"\n");
   xp=vector(1,npar);    xp=vector(1,npar);
   dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);    dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
   doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));    doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath));
      mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage);
     varpij=ma3x(1,nlstate*(nlstate+ndeath),1,nlstate*(nlstate+ndeath),(int) bage, (int) fage);
     first=1;
     if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {
       printf("Problem with gnuplot file: %s\n", optionfilegnuplot);
       exit(0);
     }
     else{
       fprintf(ficgp,"\n# Routine varprob");
     }
     if((fichtm=fopen(optionfilehtm,"a"))==NULL) {
       printf("Problem with html file: %s\n", optionfilehtm);
       exit(0);
     }
     else{
       fprintf(fichtm,"\n<H2> Computing matrix of variance-covariance of step probabilities</h2>\n");
       fprintf(fichtm,"\n<br> We have drawn ellipsoids of confidence around the p<inf>ij</inf>, p<inf>kl</inf> to understand the covariance between two incidences. They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");
       fprintf(fichtm,"\n<br> We have drawn x'cov<sup>-1</sup>x = 4 where x is the column vector (pij,pkl). It means that if pij and pkl where uncorrelated the (2X2) matrix would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 standard deviations wide on each axis. <br> When both incidences are correlated we diagonalised the inverse of the covariance matrix and made the appropriate rotation.<br> \n");
   
     }
   cov[1]=1;    cov[1]=1;
   j=cptcoveff;    j=cptcoveff;
   if (cptcovn<1) {j=1;ncodemax[1]=1;}    if (cptcovn<1) {j=1;ncodemax[1]=1;}
Line 2003  fprintf(ficresprob,"#One-step probabilit Line 2027  fprintf(ficresprob,"#One-step probabilit
   
     if  (cptcovn>0) {      if  (cptcovn>0) {
       fprintf(ficresprob, "\n#********** Variable ");        fprintf(ficresprob, "\n#********** Variable ");
         fprintf(ficresprobcov, "\n#********** Variable ");
         fprintf(ficgp, "\n#********** Variable ");
         fprintf(fichtm, "\n<h4>********** Variable</h4>\n ");
         fprintf(ficresprobcor, "\n#********** Variable ");
       for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);        for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
       fprintf(ficresprob, "**********\n#");        fprintf(ficresprob, "**********\n#");
         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
         fprintf(ficresprobcov, "**********\n#");
         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
         fprintf(ficgp, "**********\n#");
         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, "# V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
         fprintf(ficgp, "**********\n#");
         for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
         fprintf(fichtm, "**********\n#");
     }      }
         
       for (age=bage; age<=fage; age ++){        for (age=bage; age<=fage; age ++){
         cov[2]=age;          cov[2]=age;
         for (k=1; k<=cptcovn;k++) {          for (k=1; k<=cptcovn;k++) {
           cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];            cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];
            
         }          }
         for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];          for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
         for (k=1; k<=cptcovprod;k++)          for (k=1; k<=cptcovprod;k++)
           cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];            cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
                 
         gradg=matrix(1,npar,1,9);          gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath));
         trgradg=matrix(1,9,1,npar);          trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
         gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));          gp=vector(1,(nlstate)*(nlstate+ndeath));
         gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));          gm=vector(1,(nlstate)*(nlstate+ndeath));
         
         for(theta=1; theta <=npar; theta++){          for(theta=1; theta <=npar; theta++){
           for(i=1; i<=npar; i++)            for(i=1; i<=npar; i++)
Line 2029  fprintf(ficresprob,"#One-step probabilit Line 2064  fprintf(ficresprob,"#One-step probabilit
           pmij(pmmij,cov,ncovmodel,xp,nlstate);            pmij(pmmij,cov,ncovmodel,xp,nlstate);
                     
           k=0;            k=0;
           for(i=1; i<= (nlstate+ndeath); i++){            for(i=1; i<= (nlstate); i++){
             for(j=1; j<=(nlstate+ndeath);j++){              for(j=1; j<=(nlstate+ndeath);j++){
               k=k+1;                k=k+1;
               gp[k]=pmmij[i][j];                gp[k]=pmmij[i][j];
Line 2041  fprintf(ficresprob,"#One-step probabilit Line 2076  fprintf(ficresprob,"#One-step probabilit
         
           pmij(pmmij,cov,ncovmodel,xp,nlstate);            pmij(pmmij,cov,ncovmodel,xp,nlstate);
           k=0;            k=0;
           for(i=1; i<=(nlstate+ndeath); i++){            for(i=1; i<=(nlstate); i++){
             for(j=1; j<=(nlstate+ndeath);j++){              for(j=1; j<=(nlstate+ndeath);j++){
               k=k+1;                k=k+1;
               gm[k]=pmmij[i][j];                gm[k]=pmmij[i][j];
             }              }
           }            }
             
           for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++)            for(i=1; i<= (nlstate)*(nlstate+ndeath); i++)
             gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];                gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];  
         }          }
   
         for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)          for(j=1; j<=(nlstate)*(nlstate+ndeath);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];
                 
         matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);          matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov);
         matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);          matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg);
                 
         pmij(pmmij,cov,ncovmodel,x,nlstate);          pmij(pmmij,cov,ncovmodel,x,nlstate);
                 
         k=0;          k=0;
         for(i=1; i<=(nlstate+ndeath); i++){          for(i=1; i<=(nlstate); i++){
           for(j=1; j<=(nlstate+ndeath);j++){            for(j=1; j<=(nlstate+ndeath);j++){
             k=k+1;              k=k+1;
             gm[k]=pmmij[i][j];              mu[k][(int) age]=pmmij[i][j];
           }            }
         }          }
                for(i=1;i<=(nlstate)*(nlstate+ndeath);i++)
      /*printf("\n%d ",(int)age);            for(j=1;j<=(nlstate)*(nlstate+ndeath);j++)
      for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){              varpij[i][j][(int)age] = doldm[i][j];
   
           /*printf("\n%d ",(int)age);
        for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
        printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));         printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
      }*/       }*/
   
         fprintf(ficresprob,"\n%d ",(int)age);          fprintf(ficresprob,"\n%d ",(int)age);
           fprintf(ficresprobcov,"\n%d ",(int)age);
           fprintf(ficresprobcor,"\n%d ",(int)age);
   
         for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)          for (i=1; i<=(nlstate)*(nlstate+ndeath);i++)
           fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i]));            fprintf(ficresprob,"%11.3e (%11.3e) ",mu[i][(int) age],sqrt(varpij[i][i][(int)age]));
            for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
       }            fprintf(ficresprobcov,"%11.3e ",mu[i][(int) age]);
     }            fprintf(ficresprobcor,"%11.3e ",mu[i][(int) age]);
           }
           i=0;
           for (k=1; k<=(nlstate);k++){
             for (l=1; l<=(nlstate+ndeath);l++){
               i=i++;
               fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l);
               fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l);
               for (j=1; j<=i;j++){
                 fprintf(ficresprobcov," %11.3e",varpij[i][j][(int)age]);
                 fprintf(ficresprobcor," %11.3e",varpij[i][j][(int) age]/sqrt(varpij[i][i][(int) age])/sqrt(varpij[j][j][(int)age]));
               }
             }
           }/* end of loop for state */
         } /* end of loop for age */
           /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/
         for (k1=1; k1<=(nlstate);k1++){
           for (l1=1; l1<=(nlstate+ndeath);l1++){
             if(l1==k1) continue;
             i=(k1-1)*(nlstate+ndeath)+l1;
             for (k2=1; k2<=(nlstate);k2++){
               for (l2=1; l2<=(nlstate+ndeath);l2++){
                 if(l2==k2) continue;
                 j=(k2-1)*(nlstate+ndeath)+l2;
                 if(j<=i) continue;
                 for (age=bage; age<=fage; age ++){
                   if ((int)age %5==0){
                     v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM;
                     v2=varpij[j][j][(int)age]/stepm*YEARM/stepm*YEARM;
                     cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM;
                     mu1=mu[i][(int) age]/stepm*YEARM ;
                     mu2=mu[j][(int) age]/stepm*YEARM;
                     /* Computing eigen value of matrix of covariance */
                     lc1=(v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));
                     lc2=(v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));
                     printf("Var %.4e %.4e cov %.4e Eigen %.3e %.3e\n",v1,v2,cv12,lc1,lc2);
                     /* Eigen vectors */
                     v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
                     v21=sqrt(1.-v11*v11);
                     v12=-v21;
                     v22=v11;
                     /*printf(fignu*/
                     /* mu1+ v11*lc1*cost + v12*lc2*sin(t) */
                     /* mu2+ v21*lc1*cost + v21*lc2*sin(t) */
                     if(first==1){
                       first=0;
                       fprintf(ficgp,"\nset parametric;set nolabel");
                       fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k2,l2,k1,l1);
                       fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
                       fprintf(fichtm,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup> :<a href=\"varpijgr%s%1d%1d-%1d%1d.png\">varpijgr%s%1d%1d-%1d%1d.png</A>, ",k2,l2,k1,l1,optionfilefiname,k2,l2,k1,l1,optionfilefiname,k2,l2,k1,l1);
                       fprintf(fichtm,"\n<br><img src=\"varpijgr%s%1d%1d-%1d%1d.png\">, ",optionfilefiname,k2,l2,k1,l1);
                       fprintf(ficgp,"\nset out \"varpijgr%s%1d%1d-%1d%1d.png\"",optionfilefiname,k2,l2,k1,l1);
                       fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);
                       fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);
                       fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) t \"%d\"",\
                               mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
                               mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);
                     }else{
                       first=0;
                       fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);
                       fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);
                       fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) t \"%d\"",\
                               mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
                               mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);
                     }/* if first */
                   } /* age mod 5 */
                 } /* end loop age */
                 fprintf(ficgp,"\nset out \"varpijgr%s%1d%1d-%1d%1d.png\";replot;",optionfilefiname,k2,l2,k1,l1);
                 first=1;
               } /*l12 */
             } /* k12 */
           } /*l1 */
         }/* k1 */
       } /* loop covariates */
       free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
     free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));      free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
     free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));      free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
       free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
     free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);      free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
     free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);      free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
   }    }
   free_vector(xp,1,npar);    free_vector(xp,1,npar);
   fclose(ficresprob);    fclose(ficresprob);
      fclose(ficresprobcov);
     fclose(ficresprobcor);
     fclose(ficgp);
     fclose(fichtm);
 }  }
   
   
 /******************* Printing html file ***********/  /******************* Printing html file ***********/
 void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \  void printinghtml(char fileres[], 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 optionfile[], \                    int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
  char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\                    int popforecast, int estepm ,\
  char version[], int popforecast, int estepm ){                    double jprev1, double mprev1,double anprev1, \
                     double jprev2, double mprev2,double anprev2){
   int jj1, k1, i1, cpt;    int jj1, k1, i1, cpt;
   FILE *fichtm;  
   /*char optionfilehtm[FILENAMELENGTH];*/    /*char optionfilehtm[FILENAMELENGTH];*/
     if((fichtm=fopen(optionfilehtm,"a"))==NULL)    {
   strcpy(optionfilehtm,optionfile);  
   strcat(optionfilehtm,".htm");  
   if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {  
     printf("Problem with %s \n",optionfilehtm), exit(0);      printf("Problem with %s \n",optionfilehtm), exit(0);
   }    }
   
  fprintf(fichtm,"<body> <font size=\"2\">%s </font> <hr size=\"2\" color=\"#EC5E5E\"> \n     fprintf(fichtm,"<ul><li>Result files (first order: no variance)<br>\n
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n   - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"p%s\">p%s</a> <br>\n
 \n   - Estimated transition probabilities over %d (stepm) months: <a href=\"pij%s\">pij%s</a><br>\n
 Total number of observations=%d <br>\n   - Stable prevalence in each health state: <a href=\"pl%s\">pl%s</a> <br>\n
 Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n   - Life expectancies by age and initial health status (estepm=%2d months):
 <hr  size=\"2\" color=\"#EC5E5E\">     <a href=\"e%s\">e%s</a> <br>\n</li>", \
  <ul><li>Outputs files<br>\n    jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,fileres,fileres,stepm,fileres,fileres,fileres,fileres,estepm,fileres,fileres);
  - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n  
  - Gnuplot file name: <a href=\"%s\">%s</a><br>\n   fprintf(fichtm,"\n<li> Result files (second order: variances)<br>\n
  - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n   - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n
  - Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>\n   - Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n
  - Transition probabilities: <a href=\"pij%s\">pij%s</a><br>\n   - Variance-covariance of one-step probabilities: <a href=\"probcov%s\">probcov%s</a> <br>\n
  - Life expectancies by age and initial health status (estepm=%2d months): <a href=\"e%s\">e%s</a> <br>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot,fileres,fileres,fileres,fileres,fileres,fileres,estepm,fileres,fileres);   - Correlation matrix of one-step probabilities: <a href=\"probcor%s\">probcor%s</a> <br>\n
    - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n
  fprintf(fichtm,"\n   - Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n
  - Parameter file with estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>\n   - Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
   - Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n  
  - Variances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n  
  - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>\n  
  - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);  
   
  if(popforecast==1) fprintf(fichtm,"\n   if(popforecast==1) fprintf(fichtm,"\n
  - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n   - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n
Line 2142  fprintf(fichtm," <li>Graphs</li><p>"); Line 2254  fprintf(fichtm," <li>Graphs</li><p>");
  jj1=0;   jj1=0;
  for(k1=1; k1<=m;k1++){   for(k1=1; k1<=m;k1++){
    for(i1=1; i1<=ncodemax[k1];i1++){     for(i1=1; i1<=ncodemax[k1];i1++){
        jj1++;       jj1++;
        if (cptcovn > 0) {       if (cptcovn > 0) {
          fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");         fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
          for (cpt=1; cpt<=cptcoveff;cpt++)         for (cpt=1; cpt<=cptcoveff;cpt++)
            fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);           fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
          fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
        }       }
        fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>       /* Pij */
 <img src=\"pe%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);           fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png<br>
   <img src=\"pe%s%d1.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);    
        /* Quasi-incidences */
        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: pe%s%d2.png<br>
   <img src=\"pe%s%d2.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
          /* Stable prevalence in each health state */
        for(cpt=1; cpt<nlstate;cpt++){         for(cpt=1; cpt<nlstate;cpt++){
          fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.gif<br>           fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br>
 <img src=\"p%s%d%d.gif\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  <img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
        }         }
     for(cpt=1; cpt<=nlstate;cpt++) {      for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident         fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
 interval) in state (%d): v%s%d%d.gif <br>  interval) in state (%d): v%s%d%d.png <br>
 <img src=\"v%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);    <img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),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): exp%s%d%d.gif <br>          fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br>
 <img src=\"exp%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  <img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
      }       }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and       fprintf(fichtm,"\n<br>- Total life expectancy by age and
 health expectancies in states (1) and (2): e%s%d.gif<br>  health expectancies in states (1) and (2): e%s%d.png<br>
 <img src=\"e%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);  <img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
 fprintf(fichtm,"\n</body>");  
    }  
    }     }
    }
 fclose(fichtm);  fclose(fichtm);
 }  }
   
 /******************* Gnuplot file **************/  /******************* Gnuplot file **************/
 void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){  void printinggnuplot(char fileres[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
   
   int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;    int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
     int ng;
   strcpy(optionfilegnuplot,optionfilefiname);    if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {
   strcat(optionfilegnuplot,".gp.txt");  
   if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {  
     printf("Problem with file %s",optionfilegnuplot);      printf("Problem with file %s",optionfilegnuplot);
   }    }
   
Line 2193  m=pow(2,cptcoveff); Line 2307  m=pow(2,cptcoveff);
   for (cpt=1; cpt<= nlstate ; cpt ++) {    for (cpt=1; cpt<= nlstate ; cpt ++) {
    for (k1=1; k1<= m ; k1 ++) {     for (k1=1; k1<= m ; k1 ++) {
   
      fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);  #ifdef windows
        fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);
   #endif
   #ifdef unix
   fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
   fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres);
   #endif
   
 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)");
Line 2210  for (i=1; i<= nlstate ; i ++) { Line 2331  for (i=1; i<= nlstate ; i ++) {
   else fprintf(ficgp," \%%*lf (\%%*lf)");    else fprintf(ficgp," \%%*lf (\%%*lf)");
 }    }  
      fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));       fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));
   #ifdef unix
 fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);  fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\n");
   #endif
    }     }
   }    }
   /*2 eme*/    /*2 eme*/
   
   for (k1=1; k1<= m ; k1 ++) {    for (k1=1; k1<= m ; k1 ++) {
     fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",ageminpar,fage);      fprintf(ficgp,"\nset out \"e%s%d.png\" \n",strtok(optionfile, "."),k1);
       fprintf(ficgp,"set ylabel \"Years\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] ",ageminpar,fage);
         
     for (i=1; i<= nlstate+1 ; i ++) {      for (i=1; i<= nlstate+1 ; i ++) {
       k=2*i;        k=2*i;
Line 2242  fprintf(ficgp,"\nset out \"v%s%d%d.gif\" Line 2365  fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
       if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");        if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");
       else fprintf(ficgp,"\" t\"\" w l 0,");        else fprintf(ficgp,"\" t\"\" w l 0,");
     }      }
     fprintf(ficgp,"\nset out \"e%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),k1);  
   }    }
     
   /*3eme*/    /*3eme*/
Line 2250  fprintf(ficgp,"\nset out \"v%s%d%d.gif\" Line 2372  fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
   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);
       fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt);        fprintf(ficgp,"\nset out \"exp%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
         fprintf(ficgp,"set ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt);
       /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);        /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
  for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");   for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
 fprintf(ficgp,"\" t \"e%d1\" w l",cpt);  fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
Line 2263  fprintf(ficgp,"\" t \"e%d1\" w l",cpt); Line 2386  fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
         fprintf(ficgp," ,\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+2*i,cpt,i+1);          fprintf(ficgp," ,\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+2*i,cpt,i+1);
   
       }        }
       fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);  
     }  
     }      }
     }
     
   /* CV preval stat */    /* CV preval stat */
     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=3;        k=3;
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);        fprintf(ficgp,"\nset out \"p%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);
   
       for (i=1; i< nlstate ; i ++)        for (i=1; i< nlstate ; i ++)
         fprintf(ficgp,"+$%d",k+i+1);          fprintf(ficgp,"+$%d",k+i+1);
Line 2284  fprintf(ficgp,"\" t \"e%d1\" w l",cpt); Line 2407  fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
         fprintf(ficgp,"+$%d",l+i+1);          fprintf(ficgp,"+$%d",l+i+1);
       }        }
       fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);          fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);  
       fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);  
     }      }
   }      }  
     
Line 2300  fprintf(ficgp,"\" t \"e%d1\" w l",cpt); Line 2422  fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
         }          }
       }        }
     }      }
     }     }
   
     for(jk=1; jk <=m; jk++) {     for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/
   fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);       for(jk=1; jk <=m; jk++) {
    i=1;         fprintf(ficgp,"\nset out \"pe%s%d%d.png\" \n",strtok(optionfile, "."),jk,ng);
    for(k2=1; k2<=nlstate; k2++) {         if (ng==2)
      k3=i;           fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
      for(k=1; k<=(nlstate+ndeath); k++) {         else
        if (k != k2){           fprintf(ficgp,"\nset title \"Probability\"\n");
         fprintf(ficgp," exp(p%d+p%d*x",i,i+1);         fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
 ij=1;         i=1;
         for(j=3; j <=ncovmodel; j++) {         for(k2=1; k2<=nlstate; k2++) {
           if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {           k3=i;
             fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);           for(k=1; k<=(nlstate+ndeath); k++) {
             ij++;             if (k != k2){
           }               if(ng==2)
           else                 fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
           fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);               else
         }                 fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
           fprintf(ficgp,")/(1");               ij=1;
                       for(j=3; j <=ncovmodel; j++) {
         for(k1=1; k1 <=nlstate; k1++){                   if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
           fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);                   fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
 ij=1;                   ij++;
           for(j=3; j <=ncovmodel; j++){                 }
           if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {                 else
             fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);                   fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
             ij++;               }
           }               fprintf(ficgp,")/(1");
           else               
             fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);               for(k1=1; k1 <=nlstate; k1++){  
           }                 fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
           fprintf(ficgp,")");                 ij=1;
         }                 for(j=3; j <=ncovmodel; j++){
         fprintf(ficgp,") t \"p%d%d\" ", k2,k);                   if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
         if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");                     fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
         i=i+ncovmodel;                     ij++;
                    }
                    else
                      fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
                  }
                  fprintf(ficgp,")");
                }
                fprintf(ficgp,") t \"p%d%d\" ", k2,k);
                if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
                i=i+ncovmodel;
              }
            }
        }         }
      }       }
    }     }
    fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);     fclose(ficgp);
    }  
      
   fclose(ficgp);  
 }  /* end gnuplot */  }  /* end gnuplot */
   
   
Line 2643  int main(int argc, char *argv[]) Line 2773  int main(int argc, char *argv[])
   double ***p3mat;    double ***p3mat;
   int *indx;    int *indx;
   char line[MAXLINE], linepar[MAXLINE];    char line[MAXLINE], linepar[MAXLINE];
   char title[MAXLINE];  
   char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];  
   char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];  
    
   char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];  
   
   char filerest[FILENAMELENGTH];  
   char fileregp[FILENAMELENGTH];  
   char popfile[FILENAMELENGTH];  
   char path[80],pathc[80],pathcd[80],pathtot[80],model[20];    char path[80],pathc[80],pathcd[80],pathtot[80],model[20];
   int firstobs=1, lastobs=10;    int firstobs=1, lastobs=10;
   int sdeb, sfin; /* Status at beginning and end */    int sdeb, sfin; /* Status at beginning and end */
Line 2679  int main(int argc, char *argv[]) Line 2800  int main(int argc, char *argv[])
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;    double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
     
   
   char version[80]="Imach version 0.8a1, June 2003, INED-EUROREVES ";  
   char *alph[]={"a","a","b","c","d","e"}, str[4];    char *alph[]={"a","a","b","c","d","e"}, str[4];
   
   
Line 3274  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3394  while((c=getc(ficpar))=='#' && c!= EOF){
  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
   
 /*------------ gnuplot -------------*/  /*------------ gnuplot -------------*/
  printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, ageminpar,agemaxpar,fage, pathc,p);    strcpy(optionfilegnuplot,optionfilefiname);
     strcat(optionfilegnuplot,".gp");
     if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
       printf("Problem with file %s",optionfilegnuplot);
     }
     fclose(ficgp);
    printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p);
   /*--------- index.htm --------*/
   
     strcpy(optionfilehtm,optionfile);
     strcat(optionfilehtm,".htm");
     if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {
       printf("Problem with %s \n",optionfilehtm), exit(0);
     }
   
     fprintf(fichtm,"<body> <font size=\"2\">%s </font> <hr size=\"2\" color=\"#EC5E5E\"> \n
   Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n
   \n
   Total number of observations=%d <br>\n
   Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
   <hr  size=\"2\" color=\"#EC5E5E\">
    <ul><li>Parameter files<br>\n
    - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n
    - Gnuplot file name: <a href=\"%s\">%s</a><br></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot);
     fclose(fichtm);
   
    printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
     
 /*------------ free_vector  -------------*/  /*------------ free_vector  -------------*/
  chdir(path);   chdir(path);
Line 3288  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3434  while((c=getc(ficpar))=='#' && c!= EOF){
  fclose(ficparo);   fclose(ficparo);
  fclose(ficres);   fclose(ficres);
   
 /*--------- index.htm --------*/  
   
   printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast,estepm);  
   
    
   /*--------------- Prevalence limit --------------*/    /*--------------- Prevalence limit --------------*/
     
   strcpy(filerespl,"pl");    strcpy(filerespl,"pl");
Line 3387  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3529  while((c=getc(ficpar))=='#' && c!= EOF){
     }      }
   }    }
   
   varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax);    varprob(optionfilefiname, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax);
   
   fclose(ficrespij);    fclose(ficrespij);
   
Line 3543  free_matrix(mint,1,maxwav,1,n); Line 3685  free_matrix(mint,1,maxwav,1,n);
   free_matrix(agev,1,maxwav,1,imx);    free_matrix(agev,1,maxwav,1,imx);
   free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);    free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
   
     fprintf(fichtm,"\n</body>");
     fclose(fichtm);
     fclose(ficgp);
    
   
   if(erreur >0)    if(erreur >0)
     printf("End of Imach with error or warning %d\n",erreur);      printf("End of Imach with error or warning %d\n",erreur);
   else   printf("End of Imach\n");    else   printf("End of Imach\n");
Line 3554  free_matrix(mint,1,maxwav,1,n); Line 3701  free_matrix(mint,1,maxwav,1,n);
   
   
  end:   end:
   #ifdef windows
   /* chdir(pathcd);*/    /* chdir(pathcd);*/
   #endif
  /*system("wgnuplot graph.plt");*/   /*system("wgnuplot graph.plt");*/
  /*system("../gp37mgw/wgnuplot graph.plt");*/   /*system("../gp37mgw/wgnuplot graph.plt");*/
  /*system("cd ../gp37mgw");*/   /*system("cd ../gp37mgw");*/
Line 3564  free_matrix(mint,1,maxwav,1,n); Line 3713  free_matrix(mint,1,maxwav,1,n);
  strcat(plotcmd,optionfilegnuplot);   strcat(plotcmd,optionfilegnuplot);
  system(plotcmd);   system(plotcmd);
   
  /*#ifdef windows*/  #ifdef windows
   while (z[0] != 'q') {    while (z[0] != 'q') {
     /* chdir(path); */      /* chdir(path); */
     printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");      printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");
Line 3574  free_matrix(mint,1,maxwav,1,n); Line 3723  free_matrix(mint,1,maxwav,1,n);
     else if (z[0] == 'g') system(plotcmd);      else if (z[0] == 'g') system(plotcmd);
     else if (z[0] == 'q') exit(0);      else if (z[0] == 'q') exit(0);
   }    }
   /*#endif */  #endif
 }  }
   
   
Removed from v.1.41.2.2  
changed lines
  Added in v.1.48

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