Diff for /imach/src/imach.c between versions 1.89 and 1.96

version 1.89, 2003/06/24 12:30:52 version 1.96, 2003/07/15 15:38:55
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.96  2003/07/15 15:38:55  brouard
     * imach.c (Repository): Errors in subdirf, 2, 3 while printing tmpout is
     rewritten within the same printf. Workaround: many printfs.
   
     Revision 1.95  2003/07/08 07:54:34  brouard
     * imach.c (Repository):
     (Repository): Using imachwizard code to output a more meaningful covariance
     matrix (cov(a12,c31) instead of numbers.
   
     Revision 1.94  2003/06/27 13:00:02  brouard
     Just cleaning
   
     Revision 1.93  2003/06/25 16:33:55  brouard
     (Module): On windows (cygwin) function asctime_r doesn't
     exist so I changed back to asctime which exists.
     (Module): Version 0.96b
   
     Revision 1.92  2003/06/25 16:30:45  brouard
     (Module): On windows (cygwin) function asctime_r doesn't
     exist so I changed back to asctime which exists.
   
     Revision 1.91  2003/06/25 15:30:29  brouard
     * imach.c (Repository): Duplicated warning errors corrected.
     (Repository): Elapsed time after each iteration is now output. It
     helps to forecast when convergence will be reached. Elapsed time
     is stamped in powell.  We created a new html file for the graphs
     concerning matrix of covariance. It has extension -cov.htm.
   
     Revision 1.90  2003/06/24 12:34:15  brouard
     (Module): Some bugs corrected for windows. Also, when
     mle=-1 a template is output in file "or"mypar.txt with the design
     of the covariance matrix to be input.
   
   Revision 1.89  2003/06/24 12:30:52  brouard    Revision 1.89  2003/06/24 12:30:52  brouard
   (Module): Some bugs corrected for windows. Also, when    (Module): Some bugs corrected for windows. Also, when
   mle=-1 a template is output in file "or"mypar.txt with the design    mle=-1 a template is output in file "or"mypar.txt with the design
Line 141 Line 174
 #include <time.h>  #include <time.h>
 #include "timeval.h"  #include "timeval.h"
   
   /* #include <libintl.h> */
   /* #define _(String) gettext (String) */
   
 #define MAXLINE 256  #define MAXLINE 256
 #define GNUPLOTPROGRAM "gnuplot"  #define GNUPLOTPROGRAM "gnuplot"
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/  /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
Line 172 Line 208
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
   
 char version[]="Imach version 0.96a, June 2003, INED-EUROREVES ";  char version[]="Imach version 0.96c, July 2003, INED-EUROREVES ";
 char fullversion[]="$Revision$ $Date$";   char fullversion[]="$Revision$ $Date$"; 
 int erreur; /* Error number */  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;
 int npar=NPARMAX;  int npar=NPARMAX;
Line 206  char fileresilk[FILENAMELENGTH]; /* File Line 242  char fileresilk[FILENAMELENGTH]; /* File
 FILE *ficresilk;  FILE *ficresilk;
 FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;  FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;
 FILE *ficresprobmorprev;  FILE *ficresprobmorprev;
 FILE *fichtm; /* Html File */  FILE *fichtm, *fichtmcov; /* Html File */
 FILE *ficreseij;  FILE *ficreseij;
 char filerese[FILENAMELENGTH];  char filerese[FILENAMELENGTH];
 FILE  *ficresvij;  FILE  *ficresvij;
Line 216  char fileresvpl[FILENAMELENGTH]; Line 252  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];  char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
 char tmpout[FILENAMELENGTH];   char tmpout[FILENAMELENGTH],  tmpout2[FILENAMELENGTH]; 
 char command[FILENAMELENGTH];  char command[FILENAMELENGTH];
 int  outcmd=0;  int  outcmd=0;
   
 char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];  char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
 char lfileres[FILENAMELENGTH];  
 char filelog[FILENAMELENGTH]; /* Log file */  char filelog[FILENAMELENGTH]; /* Log file */
 char filerest[FILENAMELENGTH];  char filerest[FILENAMELENGTH];
 char fileregp[FILENAMELENGTH];  char fileregp[FILENAMELENGTH];
 char popfile[FILENAMELENGTH];  char popfile[FILENAMELENGTH];
   
 char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];  char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilehtmcov[FILENAMELENGTH] ;
   
   struct timeval start_time, end_time, curr_time, last_time, forecast_time;
   struct timezone tzp;
   extern int gettimeofday();
   struct tm tmg, tm, tmf, *gmtime(), *localtime();
   long time_value;
   extern long time();
   char strcurr[80], strfor[80];
   
 #define NR_END 1  #define NR_END 1
 #define FREE_ARG char*  #define FREE_ARG char*
Line 533  void free_ma3x(double ***m, long nrl, lo Line 577  void free_ma3x(double ***m, long nrl, lo
   free((FREE_ARG)(m+nrl-NR_END));    free((FREE_ARG)(m+nrl-NR_END));
 }  }
   
   /*************** function subdirf ***********/
   char *subdirf(char fileres[])
   {
     /* Caution optionfilefiname is hidden */
     strcpy(tmpout,optionfilefiname);
     strcat(tmpout,"/"); /* Add to the right */
     strcat(tmpout,fileres);
     return tmpout;
   }
   
   /*************** function subdirf2 ***********/
   char *subdirf2(char fileres[], char *preop)
   {
     
     /* Caution optionfilefiname is hidden */
     strcpy(tmpout,optionfilefiname);
     strcat(tmpout,"/");
     strcat(tmpout,preop);
     strcat(tmpout,fileres);
     return tmpout;
   }
   
   /*************** function subdirf3 ***********/
   char *subdirf3(char fileres[], char *preop, char *preop2)
   {
     
     /* Caution optionfilefiname is hidden */
     strcpy(tmpout,optionfilefiname);
     strcat(tmpout,"/");
     strcat(tmpout,preop);
     strcat(tmpout,preop2);
     strcat(tmpout,fileres);
     return tmpout;
   }
   
 /***************** f1dim *************************/  /***************** f1dim *************************/
 extern int ncom;   extern int ncom; 
 extern double *pcom,*xicom;  extern double *pcom,*xicom;
Line 708  void linmin(double p[], double xi[], int Line 787  void linmin(double p[], double xi[], int
   free_vector(pcom,1,n);     free_vector(pcom,1,n); 
 }   } 
   
   char *asc_diff_time(long time_sec, char ascdiff[])
   {
     long sec_left, days, hours, minutes;
     days = (time_sec) / (60*60*24);
     sec_left = (time_sec) % (60*60*24);
     hours = (sec_left) / (60*60) ;
     sec_left = (sec_left) %(60*60);
     minutes = (sec_left) /60;
     sec_left = (sec_left) % (60);
     sprintf(ascdiff,"%d day(s) %d hour(s) %d minute(s) %d second(s)",days, hours, minutes, sec_left);  
     return ascdiff;
   }
   
 /*************** powell ************************/  /*************** powell ************************/
 void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret,   void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret, 
             double (*func)(double []))               double (*func)(double [])) 
Line 718  void powell(double p[], double **xi, int Line 810  void powell(double p[], double **xi, int
   double del,t,*pt,*ptt,*xit;    double del,t,*pt,*ptt,*xit;
   double fp,fptt;    double fp,fptt;
   double *xits;    double *xits;
     int niterf, itmp;
   
   pt=vector(1,n);     pt=vector(1,n); 
   ptt=vector(1,n);     ptt=vector(1,n); 
   xit=vector(1,n);     xit=vector(1,n); 
Line 728  void powell(double p[], double **xi, int Line 822  void powell(double p[], double **xi, int
     fp=(*fret);       fp=(*fret); 
     ibig=0;       ibig=0; 
     del=0.0;       del=0.0; 
     printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret);      last_time=curr_time;
     fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f",*iter,*fret);      (void) gettimeofday(&curr_time,&tzp);
     fprintf(ficrespow,"%d %.12f",*iter,*fret);      printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout);
       fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);
       fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec);
     for (i=1;i<=n;i++) {      for (i=1;i<=n;i++) {
       printf(" %d %.12f",i, p[i]);        printf(" %d %.12f",i, p[i]);
       fprintf(ficlog," %d %.12lf",i, p[i]);        fprintf(ficlog," %d %.12lf",i, p[i]);
Line 738  void powell(double p[], double **xi, int Line 834  void powell(double p[], double **xi, int
     }      }
     printf("\n");      printf("\n");
     fprintf(ficlog,"\n");      fprintf(ficlog,"\n");
     fprintf(ficrespow,"\n");      fprintf(ficrespow,"\n");fflush(ficrespow);
       if(*iter <=3){
         tm = *localtime(&curr_time.tv_sec);
         strcpy(strcurr,asctime(&tmf));
   /*       asctime_r(&tm,strcurr); */
         forecast_time=curr_time;
         itmp = strlen(strcurr);
         if(strcurr[itmp-1]=='\n')
           strcurr[itmp-1]='\0';
         printf("\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);
         fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);
         for(niterf=10;niterf<=30;niterf+=10){
           forecast_time.tv_sec=curr_time.tv_sec+(niterf-*iter)*(curr_time.tv_sec-last_time.tv_sec);
           tmf = *localtime(&forecast_time.tv_sec);
   /*      asctime_r(&tmf,strfor); */
           strcpy(strfor,asctime(&tmf));
           itmp = strlen(strfor);
           if(strfor[itmp-1]=='\n')
           strfor[itmp-1]='\0';
           printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s or\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
           fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s or\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
         }
       }
     for (i=1;i<=n;i++) {       for (i=1;i<=n;i++) { 
       for (j=1;j<=n;j++) xit[j]=xi[j][i];         for (j=1;j<=n;j++) xit[j]=xi[j][i]; 
       fptt=(*fret);         fptt=(*fret); 
Line 1151  double func( double *x) Line 1269  double func( double *x)
           oldm=newm;            oldm=newm;
         } /* end mult */          } /* end mult */
               
         /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */  
         /* But now since version 0.9 we anticipate for bias and large stepm.  
          * If stepm is larger than one month (smallest stepm) and if the exact delay   
          * (in months) between two waves is not a multiple of stepm, we rounded to   
          * the nearest (and in case of equal distance, to the lowest) interval but now  
          * we keep into memory the bias bh[mi][i] and also the previous matrix product  
          * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the  
          * probability in order to take into account the bias as a fraction of the way  
          * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies  
          * -stepm/2 to stepm/2 .  
          * For stepm=1 the results are the same as for previous versions of Imach.  
          * For stepm > 1 the results are less biased than in previous versions.   
          */  
         s1=s[mw[mi][i]][i];          s1=s[mw[mi][i]][i];
         s2=s[mw[mi+1][i]][i];          s2=s[mw[mi+1][i]][i];
         bbh=(double)bh[mi][i]/(double)stepm;           bbh=(double)bh[mi][i]/(double)stepm; 
         /* bias is positive if real duration  
          * is higher than the multiple of stepm and negative otherwise.  
          */  
         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 */
         /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/  
         /*lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-+bh)*out[s1][s2])); */ /* exponential interpolation */  
         /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/  
         /*if(lli ==000.0)*/  
         /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][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;
Line 1202  double func( double *x) Line 1299  double func( double *x)
           oldm=newm;            oldm=newm;
         } /* end mult */          } /* end mult */
               
         /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */  
         /* But now since version 0.9 we anticipate for bias and large stepm.  
          * If stepm is larger than one month (smallest stepm) and if the exact delay   
          * (in months) between two waves is not a multiple of stepm, we rounded to   
          * the nearest (and in case of equal distance, to the lowest) interval but now  
          * we keep into memory the bias bh[mi][i] and also the previous matrix product  
          * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the  
          * probability in order to take into account the bias as a fraction of the way  
          * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies  
          * -stepm/2 to stepm/2 .  
          * For stepm=1 the results are the same as for previous versions of Imach.  
          * For stepm > 1 the results are less biased than in previous versions.   
          */  
         s1=s[mw[mi][i]][i];          s1=s[mw[mi][i]][i];
         s2=s[mw[mi+1][i]][i];          s2=s[mw[mi+1][i]][i];
         bbh=(double)bh[mi][i]/(double)stepm;           bbh=(double)bh[mi][i]/(double)stepm; 
         /* bias is positive if real duration  
          * is higher than the multiple of stepm and negative otherwise.  
          */  
         /* 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]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */          lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
         /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/  
         /*if(lli ==000.0)*/  
         /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][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;
Line 1391  double funcone( double *x) Line 1468  double funcone( double *x)
   return -l;    return -l;
 }  }
   
 char *subdirf(char fileres[])  
 {  
     
   strcpy(tmpout,optionfilefiname);  
   strcat(tmpout,"/"); /* Add to the right */  
   strcat(tmpout,fileres);  
   return tmpout;  
 }  
   
 char *subdirf2(char fileres[], char *preop)  
 {  
     
   strcpy(tmpout,optionfilefiname);  
   strcat(tmpout,"/");  
   strcat(tmpout,preop);  
   strcat(tmpout,fileres);  
   return tmpout;  
 }  
 char *subdirf3(char fileres[], char *preop, char *preop2)  
 {  
     
   strcpy(tmpout,optionfilefiname);  
   strcat(tmpout,"/");  
   strcat(tmpout,preop);  
   strcat(tmpout,preop2);  
   strcat(tmpout,fileres);  
   return tmpout;  
 }  
   
   /*************** function likelione ***********/
 void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*funcone)(double []))  void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*funcone)(double []))
 {  {
   /* This routine should help understanding what is done with     /* This routine should help understanding what is done with 
Line 2031  void  concatwav(int wav[], int **dh, int Line 2081  void  concatwav(int wav[], int **dh, int
   
     wav[i]=mi;      wav[i]=mi;
     if(mi==0){      if(mi==0){
         nbwarn++;
       if(first==0){        if(first==0){
         printf("Warning! None valid information for:%ld line=%d (skipped) and may be others, see log file\n",num[i],i);          printf("Warning! None valid information for:%ld line=%d (skipped) and may be others, see log file\n",num[i],i);
         first=1;          first=1;
Line 2051  void  concatwav(int wav[], int **dh, int Line 2102  void  concatwav(int wav[], int **dh, int
             j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12);               j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); 
             if(j==0) j=1;  /* Survives at least one month after exam */              if(j==0) j=1;  /* Survives at least one month after exam */
             else if(j<0){              else if(j<0){
                 nberr++;
               printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);                printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
               j=1; /* Careful Patch */                j=1; /* Temporary Dangerous patch */
               printf("   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n  You MUST fix the contradiction between dates.\n",stepm);                printf("   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n  You MUST fix the contradiction between dates.\n",stepm);
               printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);                fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
               fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n  You MUST fix the contradiction between dates.\n",stepm);                fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n  You MUST fix the contradiction between dates.\n",stepm);
             }              }
             k=k+1;              k=k+1;
Line 2074  void  concatwav(int wav[], int **dh, int Line 2126  void  concatwav(int wav[], int **dh, int
           /*        if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */            /*        if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
           /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/            /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
           if(j<0){            if(j<0){
               nberr++;
             printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);              printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
             fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);              fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
           }            }
Line 2250  void evsij(char fileres[], double ***eij Line 2303  void evsij(char fileres[], double ***eij
   
     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 */
   
      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 2804  void varprob(char optionfilefiname[], do Line 2857  void varprob(char optionfilefiname[], do
   fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");    fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
   fprintf(fichtm,"\n");    fprintf(fichtm,"\n");
   
   fprintf(fichtm,"\n<li><h4> Computing matrix of variance-covariance of step probabilities</h4></li>\n");    fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of pairs of step probabilities (drawings)</a></h4></li>\n",optionfilehtmcov);
   fprintf(fichtm,"\nWe 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(fichtmcov,"\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\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");    file %s<br>\n",optionfilehtmcov);
     fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated\
   and drawn. It helps understanding how is the covariance between two incidences.\
    They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");
     fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \
   It can be understood this way: if pij and pkl where uncorrelated the (2x2) matrix of covariance \
   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>\
    Now, if both incidences are correlated (usual case) we diagonalised the inverse of the covariance matrix\
    and made the appropriate rotation to look at the uncorrelated principal directions.<br>\
   To be simple, these graphs help to understand the significativity of each parameter in relation to a second other one.<br> \n");
   
   cov[1]=1;    cov[1]=1;
   tj=cptcoveff;    tj=cptcoveff;
Line 2828  void varprob(char optionfilefiname[], do Line 2891  void varprob(char optionfilefiname[], do
         fprintf(ficgp, "**********\n#\n");          fprintf(ficgp, "**********\n#\n");
                   
                   
         fprintf(fichtm, "\n<hr  size=\"2\" color=\"#EC5E5E\">********** Variable ");           fprintf(fichtmcov, "\n<hr  size=\"2\" color=\"#EC5E5E\">********** Variable "); 
         for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);          for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
         fprintf(fichtm, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");          fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");
                   
         fprintf(ficresprobcor, "\n#********** Variable ");              fprintf(ficresprobcor, "\n#********** Variable ");    
         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);          for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
Line 2988  void varprob(char optionfilefiname[], do Line 3051  void varprob(char optionfilefiname[], do
                     fprintf(ficgp,"\nset parametric;unset label");                      fprintf(ficgp,"\nset parametric;unset label");
                     fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2);                      fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2);
                     fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");                      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>\                      fprintf(fichtmcov,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\
  :<a href=\"%s%d%1d%1d-%1d%1d.png\">\   :<a href=\"%s%d%1d%1d-%1d%1d.png\">\
 %s%d%1d%1d-%1d%1d.png</A>, ",k1,l1,k2,l2,\  %s%d%1d%1d-%1d%1d.png</A>, ",k1,l1,k2,l2,\
                             subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2,\                              subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2,\
                             subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);                              subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
                     fprintf(fichtm,"\n<br><img src=\"%s%d%1d%1d-%1d%1d.png\"> ",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);                      fprintf(fichtmcov,"\n<br><img src=\"%s%d%1d%1d-%1d%1d.png\"> ",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
                     fprintf(fichtm,"\n<br> Correlation at age %d (%.3f),",(int) age, c12);                      fprintf(fichtmcov,"\n<br> Correlation at age %d (%.3f),",(int) age, c12);
                     fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\"",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);                      fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\"",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
                     fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);                      fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
                     fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);                      fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
Line 3003  void varprob(char optionfilefiname[], do Line 3066  void varprob(char optionfilefiname[], do
                             mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));                              mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));
                   }else{                    }else{
                     first=0;                      first=0;
                     fprintf(fichtm," %d (%.3f),",(int) age, c12);                      fprintf(fichtmcov," %d (%.3f),",(int) age, c12);
                     fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);                      fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
                     fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);                      fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
                     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)) not",\                      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)) not",\
Line 3026  void varprob(char optionfilefiname[], do Line 3089  void varprob(char optionfilefiname[], do
   fclose(ficresprob);    fclose(ficresprob);
   fclose(ficresprobcov);    fclose(ficresprobcov);
   fclose(ficresprobcor);    fclose(ficresprobcor);
   /*  fclose(ficgp);*/    fflush(ficgp);
     fflush(fichtmcov);
 }  }
   
   
Line 3038  void printinghtml(char fileres[], char t Line 3102  void printinghtml(char fileres[], char t
                   double jprev1, double mprev1,double anprev1, \                    double jprev1, double mprev1,double anprev1, \
                   double jprev2, double mprev2,double anprev2){                    double jprev2, double mprev2,double anprev2){
   int jj1, k1, i1, cpt;    int jj1, k1, i1, cpt;
   /*char optionfilehtm[FILENAMELENGTH];*/  
 /*   if((fichtm=fopen(optionfilehtm,"a"))==NULL)    { */  
 /*     printf("Problem with %s \n",optionfilehtm), exit(0); */  
 /*     fprintf(ficlog,"Problem with %s \n",optionfilehtm), exit(0); */  
 /*   } */  
   
    fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n \     fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n \
  - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n \   - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ",
  - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n \             jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p"));
  - Stable prevalence in each health state: <a href=\"%s\">%s</a> <br>\n \     fprintf(fichtm,"\
    - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
              stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij"));
      fprintf(fichtm,"\
    - Stable prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
              subdirf2(fileres,"pl"),subdirf2(fileres,"pl"));
      fprintf(fichtm,"\
  - Life expectancies by age and initial health status (estepm=%2d months): \   - Life expectancies by age and initial health status (estepm=%2d months): \
    <a href=\"%s\">%s</a> <br>\n</li>", \     <a href=\"%s\">%s</a> <br>\n</li>",
            jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p"),\  
            stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij"),\  
            subdirf2(fileres,"pl"),subdirf2(fileres,"pl"),\  
            estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));             estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));
   
 fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");  fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
Line 3086  fprintf(fichtm," \n<ul><li><b>Graphs</b> Line 3148  fprintf(fichtm," \n<ul><li><b>Graphs</b>
         fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): %s%d%d.png <br> \          fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): %s%d%d.png <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);
      }       }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and \  
 health expectancies in states (1) and (2): %s%d.png<br>\  
 <img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);  
    } /* end i1 */     } /* end i1 */
  }/* End k1 */   }/* End k1 */
  fprintf(fichtm,"</ul>");   fprintf(fichtm,"</ul>");
   
   
  fprintf(fichtm,"\n<br><li><h4> Result files (second order: variances)</h4>\n\   fprintf(fichtm,"\
  - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n\  \n<br><li><h4> Result files (second order: variances)</h4>\n\
  - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n\   - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n", rfileres,rfileres);
  - Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n\  
  - Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n\   fprintf(fichtm," - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
  - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"%s\">%s</a><br>\n\           subdirf2(fileres,"prob"),subdirf2(fileres,"prob"));
  - Health expectancies with their variances (no covariance): <a href=\"%s\">%s</a> <br>\n\   fprintf(fichtm,"\
    - Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
            subdirf2(fileres,"probcov"),subdirf2(fileres,"probcov"));
   
    fprintf(fichtm,"\
    - Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
            subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor"));
    fprintf(fichtm,"\
    - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"%s\">%s</a><br>\n",
            estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"));
    fprintf(fichtm,"\
    - Health expectancies with their variances (no covariance): <a href=\"%s\">%s</a> <br>\n",
            subdirf2(fileres,"t"),subdirf2(fileres,"t"));
    fprintf(fichtm,"\
  - Standard deviation of stable prevalences: <a href=\"%s\">%s</a> <br>\n",\   - Standard deviation of stable prevalences: <a href=\"%s\">%s</a> <br>\n",\
          rfileres,rfileres,\  
          subdirf2(fileres,"prob"),subdirf2(fileres,"prob"),\  
          subdirf2(fileres,"probcov"),subdirf2(fileres,"probcov"),\  
          subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor"),\  
          estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"),\  
          subdirf2(fileres,"t"),subdirf2(fileres,"t"),\  
          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 3116  health expectancies in states (1) and (2 Line 3182  health expectancies in states (1) and (2
 /*      <br>",fileres,fileres,fileres,fileres); */  /*      <br>",fileres,fileres,fileres,fileres); */
 /*  else  */  /*  else  */
 /*    fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */  /*    fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */
 fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");   fflush(fichtm);
    fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");
   
  m=cptcoveff;   m=cptcoveff;
  if (cptcovn < 1) {m=1;ncodemax[1]=1;}   if (cptcovn < 1) {m=1;ncodemax[1]=1;}
Line 3132  fprintf(fichtm," <ul><li><b>Graphs</b></ Line 3199  fprintf(fichtm," <ul><li><b>Graphs</b></
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }       }
      for(cpt=1; cpt<=nlstate;cpt++) {       for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"<br>- Observed and period prevalence (with confident\         fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \
 interval) in state (%d): %s%d%d.png <br>\  prevalence (with 95%% confidence interval) in state (%d): %s%d%d.png <br>\
 <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"vr"),cpt,jj1,subdirf2(optionfilefiname,"vr"),cpt,jj1);    <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);  
      }       }
        fprintf(fichtm,"\n<br>- Total life expectancy by age and \
   health expectancies in states (1) and (2): %s%d.png<br>\
   <img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
    } /* end i1 */     } /* end i1 */
  }/* End k1 */   }/* End k1 */
  fprintf(fichtm,"</ul>");   fprintf(fichtm,"</ul>");
Line 3669  int fileappend(FILE *fichier, char *opti Line 3739  int fileappend(FILE *fichier, char *opti
   fflush(fichier);    fflush(fichier);
   return (1);    return (1);
 }  }
   
   
   /**************** function prwizard **********************/
 void prwizard(int ncovmodel, int nlstate, int ndeath,  char model[], FILE *ficparo)  void prwizard(int ncovmodel, int nlstate, int ndeath,  char model[], FILE *ficparo)
 {  {
   
     /* Wizard to print covariance matrix template */
   
   char ca[32], cb[32], cc[32];    char ca[32], cb[32], cc[32];
   int i,j, k, l, li, lj, lk, ll, jj, npar, itimes;    int i,j, k, l, li, lj, lk, ll, jj, npar, itimes;
   int numlinepar;    int numlinepar;
Line 3785  void prwizard(int ncovmodel, int nlstate Line 3860  void prwizard(int ncovmodel, int nlstate
         } /* end k*/          } /* end k*/
       } /*end j */        } /*end j */
     } /* end i */      } /* end i */
   }    } /* end itimes */
   
 } /* end of prwizard */  } /* end of prwizard */
   
   
 /***********************************************/  /***********************************************/
 /**************** Main Program *****************/  /**************** Main Program *****************/
 /***********************************************/  /***********************************************/
Line 3797  int main(int argc, char *argv[]) Line 3873  int main(int argc, char *argv[])
 {  {
   int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);    int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
   int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;    int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;
   int jj, imk;    int jj, ll, li, lj, lk, imk;
   int numlinepar=0; /* Current linenumber of parameter file */    int numlinepar=0; /* Current linenumber of parameter file */
     int itimes;
   
     char ca[32], cb[32], cc[32];
   /*  FILE *fichtm; *//* Html File */    /*  FILE *fichtm; *//* Html File */
   /* FILE *ficgp;*/ /*Gnuplot File */    /* FILE *ficgp;*/ /*Gnuplot File */
   double agedeb, agefin,hf;    double agedeb, agefin,hf;
Line 3852  int main(int argc, char *argv[]) Line 3931  int main(int argc, char *argv[])
   int lstra;    int lstra;
   
   long total_usecs;    long total_usecs;
   struct timeval start_time, end_time, curr_time;  
   struct timezone tzp;  
   extern int gettimeofday();  
   struct tm tmg, tm, *gmtime(), *localtime();  
   long time_value;  
   extern long time();  
     
   /*   setlocale (LC_ALL, ""); */
   /*   bindtextdomain (PACKAGE, LOCALEDIR); */
   /*   textdomain (PACKAGE); */
   /*   setlocale (LC_CTYPE, ""); */
   /*   setlocale (LC_MESSAGES, ""); */
   
   /*   gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */    /*   gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
   (void) gettimeofday(&start_time,&tzp);    (void) gettimeofday(&start_time,&tzp);
     curr_time=start_time;
   tm = *localtime(&start_time.tv_sec);    tm = *localtime(&start_time.tv_sec);
   tmg = *gmtime(&start_time.tv_sec);    tmg = *gmtime(&start_time.tv_sec);
   strcpy(strstart,asctime(&tm));    strcpy(strstart,asctime(&tm));
Line 3881  int main(int argc, char *argv[]) Line 3961  int main(int argc, char *argv[])
 *  printf("tim_value=%d,asctime=%s\n",time_value,strstart);   *  printf("tim_value=%d,asctime=%s\n",time_value,strstart); 
 */  */
   
     nberr=0; /* Number of errors and warnings */
     nbwarn=0;
   getcwd(pathcd, size);    getcwd(pathcd, size);
   
   printf("\n%s\n%s",version,fullversion);    printf("\n%s\n%s",version,fullversion);
Line 3929  int main(int argc, char *argv[]) Line 4011  int main(int argc, char *argv[])
  optionfilext=%s\n\   optionfilext=%s\n\
  optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);   optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
   
   printf("Localtime (at start):%s",strstart);    printf("Local time (at start):%s",strstart);
   fprintf(ficlog,"Localtime (at start): %s",strstart);    fprintf(ficlog,"Local time (at start): %s",strstart);
   fflush(ficlog);    fflush(ficlog);
   /*   (void) gettimeofday(&curr_time,&tzp); */
   /*   printf("Elapsed time %d\n", asc_diff_time(curr_time.tv_sec-start_time.tv_sec,tmpout)); */
   
   /* */    /* */
   strcpy(fileres,"r");    strcpy(fileres,"r");
Line 4323  int main(int argc, char *argv[]) Line 4407  int main(int argc, char *argv[])
         s[m][i]=-1;          s[m][i]=-1;
       }        }
       if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){        if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){
           nberr++;
         printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);          printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);
         fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);          fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);
         s[m][i]=-1;          s[m][i]=-1;
       }        }
       if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){        if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
           nberr++;
         printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]);           printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); 
         fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]);           fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); 
         s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */          s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */
Line 4346  int main(int argc, char *argv[]) Line 4432  int main(int argc, char *argv[])
           /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/            /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
             else {              else {
               if ((int)andc[i]!=9999){                if ((int)andc[i]!=9999){
                   nbwarn++;
                 printf("Warning negative age at death: %ld line:%d\n",num[i],i);                  printf("Warning negative age at death: %ld line:%d\n",num[i],i);
                 fprintf(ficlog,"Warning negative age at death: %ld line:%d\n",num[i],i);                  fprintf(ficlog,"Warning negative age at death: %ld line:%d\n",num[i],i);
                 agev[m][i]=-1;                  agev[m][i]=-1;
Line 4382  int main(int argc, char *argv[]) Line 4469  int main(int argc, char *argv[])
   for (i=1; i<=imx; i++)  {    for (i=1; i<=imx; i++)  {
     for(m=firstpass; (m<=lastpass); m++){      for(m=firstpass; (m<=lastpass); m++){
       if (s[m][i] > (nlstate+ndeath)) {        if (s[m][i] > (nlstate+ndeath)) {
           nberr++;
         printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);               printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);     
         fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);               fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);     
         goto end;          goto end;
Line 4464  int main(int argc, char *argv[]) Line 4552  int main(int argc, char *argv[])
   /*  fclose(ficgp);*/    /*  fclose(ficgp);*/
   /*--------- index.htm --------*/    /*--------- index.htm --------*/
   
   strcpy(optionfilehtm,optionfilefiname);    strcpy(optionfilehtm,optionfilefiname); /* Main html file */
   strcat(optionfilehtm,".htm");    strcat(optionfilehtm,".htm");
   if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {    if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {
     printf("Problem with %s \n",optionfilehtm), exit(0);      printf("Problem with %s \n",optionfilehtm), exit(0);
   }    }
   
     strcpy(optionfilehtmcov,optionfilefiname); /* Only for matrix of covariance */
     strcat(optionfilehtmcov,"-cov.htm");
     if((fichtmcov=fopen(optionfilehtmcov,"w"))==NULL)    {
       printf("Problem with %s \n",optionfilehtmcov), exit(0);
     }
     else{
     fprintf(fichtmcov,"<body>\n<title>IMaCh Cov %s</title>\n <font size=\"2\">%s <br> %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",\
             fileres,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,"<body>\n<title>IMaCh %s</title>\n <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\
Line 4480  Title=%s <br>Datafile=%s Firstpass=%d La Line 4580  Title=%s <br>Datafile=%s Firstpass=%d La
  - 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,\            fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\
           model,fileres,fileres,\            fileres,fileres,\
           filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);            filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);
   /*fclose(fichtm);*/  
   fflush(fichtm);    fflush(fichtm);
   
   strcpy(pathr,path);    strcpy(pathr,path);
   strcat(pathr,optionfilefiname);    strcat(pathr,optionfilefiname);
   chdir(optionfilefiname); /* Move to directory named optionfile */    chdir(optionfilefiname); /* Move to directory named optionfile */
   strcpy(lfileres,fileres);  
   strcat(lfileres,"/");  
   strcat(lfileres,optionfilefiname);  
       
   /* Calculates basic frequencies. Computes observed prevalence at single age    /* Calculates basic frequencies. Computes observed prevalence at single age
      and prints on file fileres'p'. */       and prints on file fileres'p'. */
Line 4533  Interval (in months) between two waves: Line 4629  Interval (in months) between two waves:
   fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);    fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
       
   
   jk=1;  
   fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");    fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
   printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");    printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
   fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");    fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
   for(i=1,jk=1; i <=nlstate; i++){    for(i=1,jk=1; i <=nlstate; i++){
     for(k=1; k <=(nlstate+ndeath); k++){      for(k=1; k <=(nlstate+ndeath); k++){
       if (k != i)         if (k != i) {
         {          printf("%d%d ",i,k);
           printf("%d%d ",i,k);          fprintf(ficlog,"%d%d ",i,k);
           fprintf(ficlog,"%d%d ",i,k);          fprintf(ficres,"%1d%1d ",i,k);
           fprintf(ficres,"%1d%1d ",i,k);          for(j=1; j <=ncovmodel; j++){
           for(j=1; j <=ncovmodel; j++){            printf("%f ",p[jk]);
             printf("%f ",p[jk]);            fprintf(ficlog,"%f ",p[jk]);
             fprintf(ficlog,"%f ",p[jk]);            fprintf(ficres,"%f ",p[jk]);
             fprintf(ficres,"%f ",p[jk]);            jk++; 
             jk++;   
           }  
           printf("\n");  
           fprintf(ficlog,"\n");  
           fprintf(ficres,"\n");  
         }          }
           printf("\n");
           fprintf(ficlog,"\n");
           fprintf(ficres,"\n");
         }
     }      }
   }    }
   if(mle!=0){    if(mle!=0){
Line 4584  Interval (in months) between two waves: Line 4678  Interval (in months) between two waves:
   }    }
         
   fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");    fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
   if(mle==1)    if(mle>=1)
     printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");      printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
   fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");    fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
   for(i=1,k=1;i<=npar;i++){  /* # 121 Var(a12)\n\ */
     /*  if (k>nlstate) k=1;  /* # 122 Cov(b12,a12) Var(b12)\n\ */
         i1=(i-1)/(ncovmodel*nlstate)+1;   /* # 131 Cov(a13,a12) Cov(a13,b12, Var(a13)\n\ */
         fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);  /* # 132 Cov(b13,a12) Cov(b13,b12, Cov(b13,a13) Var(b13)\n\ */
         printf("%s%d%d",alph[k],i1,tab[i]);  /* # 212 Cov(a21,a12) Cov(a21,b12, Cov(a21,a13) Cov(a21,b13) Var(a21)\n\ */
     */  /* # 212 Cov(b21,a12) Cov(b21,b12, Cov(b21,a13) Cov(b21,b13) Cov(b21,a21) Var(b21)\n\ */
     fprintf(ficres,"%3d",i);  /* # 232 Cov(a23,a12) Cov(a23,b12, Cov(a23,a13) Cov(a23,b13) Cov(a23,a21) Cov(a23,b21) Var(a23)\n\ */
     if(mle==1)  /* # 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n" */
       printf("%3d",i);  
     fprintf(ficlog,"%3d",i);  
     for(j=1; j<=i;j++){  /* Just to have a covariance matrix which will be more understandable
       fprintf(ficres," %.5e",matcov[i][j]);     even is we still don't want to manage dictionary of variables
       if(mle==1)  */
         printf(" %.5e",matcov[i][j]);    for(itimes=1;itimes<=2;itimes++){
       fprintf(ficlog," %.5e",matcov[i][j]);      jj=0;
     }      for(i=1; i <=nlstate; i++){
     fprintf(ficres,"\n");        for(j=1; j <=nlstate+ndeath; j++){
     if(mle==1)          if(j==i) continue;
       printf("\n");          for(k=1; k<=ncovmodel;k++){
     fprintf(ficlog,"\n");            jj++;
     k++;            ca[0]= k+'a'-1;ca[1]='\0';
   }            if(itimes==1){
                  if(mle>=1)
                 printf("#%1d%1d%d",i,j,k);
               fprintf(ficlog,"#%1d%1d%d",i,j,k);
               fprintf(ficres,"#%1d%1d%d",i,j,k);
             }else{
               if(mle>=1)
                 printf("%1d%1d%d",i,j,k);
               fprintf(ficlog,"%1d%1d%d",i,j,k);
               fprintf(ficres,"%1d%1d%d",i,j,k);
             }
             ll=0;
             for(li=1;li <=nlstate; li++){
               for(lj=1;lj <=nlstate+ndeath; lj++){
                 if(lj==li) continue;
                 for(lk=1;lk<=ncovmodel;lk++){
                   ll++;
                   if(ll<=jj){
                     cb[0]= lk +'a'-1;cb[1]='\0';
                     if(ll<jj){
                       if(itimes==1){
                         if(mle>=1)
                           printf(" Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj);
                         fprintf(ficlog," Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj);
                         fprintf(ficres," Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj);
                       }else{
                         if(mle>=1)
                           printf(" %.5e",matcov[jj][ll]); 
                         fprintf(ficlog," %.5e",matcov[jj][ll]); 
                         fprintf(ficres," %.5e",matcov[jj][ll]); 
                       }
                     }else{
                       if(itimes==1){
                         if(mle>=1)
                           printf(" Var(%s%1d%1d)",ca,i,j);
                         fprintf(ficlog," Var(%s%1d%1d)",ca,i,j);
                         fprintf(ficres," Var(%s%1d%1d)",ca,i,j);
                       }else{
                         if(mle>=1)
                           printf(" %.5e",matcov[jj][ll]); 
                         fprintf(ficlog," %.5e",matcov[jj][ll]); 
                         fprintf(ficres," %.5e",matcov[jj][ll]); 
                       }
                     }
                   }
                 } /* end lk */
               } /* end lj */
             } /* end li */
             if(mle>=1)
               printf("\n");
             fprintf(ficlog,"\n");
             fprintf(ficres,"\n");
             numlinepar++;
           } /* end k*/
         } /*end j */
       } /* end i */
     } /* end itimes */
   
     fflush(ficlog);
     fflush(ficres);
   
   while((c=getc(ficpar))=='#' && c!= EOF){    while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      fgets(line, MAXLINE, ficpar);
Line 4824  Interval (in months) between two waves: Line 4977  Interval (in months) between two waves:
   fclose(ficrespij);    fclose(ficrespij);
   
   probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);    probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
     for(i=1;i<=AGESUP;i++)
       for(j=1;j<=NCOVMAX;j++)
         for(k=1;k<=NCOVMAX;k++)
           probs[i][j][k]=0.;
   
   /*---------- Forecasting ------------------*/    /*---------- Forecasting ------------------*/
   /*if((stepm == 1) && (strcmp(model,".")==0)){*/    /*if((stepm == 1) && (strcmp(model,".")==0)){*/
Line 5019  ageminpar, agemax, s[lastpass][imx], age Line 5176  ageminpar, agemax, s[lastpass][imx], age
   fflush(ficgp);    fflush(ficgp);
       
   
   if(erreur >0){    if((nberr >0) || (nbwarn>0)){
     printf("End of Imach with error or warning %d\n",erreur);      printf("End of Imach with %d errors and/or %d warnings\n",nberr,nbwarn);
     fprintf(ficlog,"End of Imach with error or warning %d\n",erreur);      fprintf(ficlog,"End of Imach with %d errors and/or warnings %d\n",nberr,nbwarn);
   }else{    }else{
    printf("End of Imach\n");      printf("End of Imach\n");
    fprintf(ficlog,"End of Imach\n");      fprintf(ficlog,"End of Imach\n");
   }    }
   printf("See log file on %s\n",filelog);    printf("See log file on %s\n",filelog);
   fclose(ficlog);  
   /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */    /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */
   (void) gettimeofday(&end_time,&tzp);    (void) gettimeofday(&end_time,&tzp);
   tm = *localtime(&end_time.tv_sec);    tm = *localtime(&end_time.tv_sec);
   tmg = *gmtime(&end_time.tv_sec);    tmg = *gmtime(&end_time.tv_sec);
   strcpy(strtend,asctime(&tm));    strcpy(strtend,asctime(&tm));
   printf("Localtime at start %s\nLocaltime at end   %s",strstart, strtend);     printf("Local time at start %s\nLocaltime at end   %s",strstart, strtend); 
   fprintf(ficlog,"Localtime at start %s\nLocal time at end   %s",strstart, strtend);     fprintf(ficlog,"Local time at start %s\nLocal time at end   %s\n",strstart, strtend); 
   /*  printf("Total time used %d Sec\n", asc_time(end_time.tv_sec -start_time.tv_sec);*/    printf("Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));
   
   printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);    printf("Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec);
   fprintf(ficlog,"Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);    fprintf(ficlog,"Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));
     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>",strstart, strtend);
   fclose(fichtm);    fclose(fichtm);
     fclose(fichtmcov);
   fclose(ficgp);    fclose(ficgp);
     fclose(ficlog);
   /*------ End -----------*/    /*------ End -----------*/
   
   end:  
 #ifdef windows  
   /* chdir(pathcd);*/  
 #endif   
   chdir(path);    chdir(path);
  /*system("wgnuplot graph.plt");*/  
  /*system("../gp37mgw/wgnuplot graph.plt");*/  
  /*system("cd ../gp37mgw");*/  
  /* system("..\\gp37mgw\\wgnuplot graph.plt");*/  
   strcpy(plotcmd,GNUPLOTPROGRAM);    strcpy(plotcmd,GNUPLOTPROGRAM);
   strcat(plotcmd," ");    strcat(plotcmd," ");
   strcat(plotcmd,optionfilegnuplot);    strcat(plotcmd,optionfilegnuplot);
   printf("Starting graphs with: %s",plotcmd);fflush(stdout);    printf("Starting graphs with: %s",plotcmd);fflush(stdout);
   system(plotcmd);    if((outcmd=system(plotcmd)) != 0){
       printf(" Problem with gnuplot\n");
     }
   printf(" Wait...");    printf(" Wait...");
   
  /*#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 and q for exiting: ");
     scanf("%s",z);      scanf("%s",z);
     if (z[0] == 'c') system("./imach");  /*     if (z[0] == 'c') system("./imach"); */
     else if (z[0] == 'e') system(optionfilehtm);      if (z[0] == 'e') system(optionfilehtm);
     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 */    end:
     while (z[0] != 'q') {
       printf("\nType  q for exiting: ");
       scanf("%s",z);
     }
 }  }
   
   

Removed from v.1.89  
changed lines
  Added in v.1.96


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