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

Diff for /imach/src/imach.c between versions 1.95 and 1.102

version 1.95, 2003/07/08 07:54:34 version 1.102, 2004/09/15 17:31:30
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.102  2004/09/15 17:31:30  brouard
     Add the possibility to read data file including tab characters.
   
     Revision 1.101  2004/09/15 10:38:38  brouard
     Fix on curr_time
   
     Revision 1.100  2004/07/12 18:29:06  brouard
     Add version for Mac OS X. Just define UNIX in Makefile
   
     Revision 1.99  2004/06/05 08:57:40  brouard
     *** empty log message ***
   
     Revision 1.98  2004/05/16 15:05:56  brouard
     New version 0.97 . First attempt to estimate force of mortality
     directly from the data i.e. without the need of knowing the health
     state at each age, but using a Gompertz model: log u =a + b*age .
     This is the basic analysis of mortality and should be done before any
     other analysis, in order to test if the mortality estimated from the
     cross-longitudinal survey is different from the mortality estimated
     from other sources like vital statistic data.
   
     The same imach parameter file can be used but the option for mle should be -3.
   
     Agnès, who wrote this part of the code, tried to keep most of the
     former routines in order to include the new code within the former code.
   
     The output is very simple: only an estimate of the intercept and of
     the slope with 95% confident intervals.
   
     Current limitations:
     A) Even if you enter covariates, i.e. with the
     model= V1+V2 equation for example, the programm does only estimate a unique global model without covariates.
     B) There is no computation of Life Expectancy nor Life Table.
   
     Revision 1.97  2004/02/20 13:25:42  lievre
     Version 0.96d. Population forecasting command line is (temporarily)
     suppressed.
   
     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    Revision 1.95  2003/07/08 07:54:34  brouard
   * imach.c (Repository):    * imach.c (Repository):
   (Repository): Using imachwizard code to output a more meaningful covariance    (Repository): Using imachwizard code to output a more meaningful covariance
Line 166 Line 208
 #include <stdlib.h>  #include <stdlib.h>
 #include <unistd.h>  #include <unistd.h>
   
 #include <sys/time.h>  /* #include <sys/time.h> */
 #include <time.h>  #include <time.h>
 #include "timeval.h"  #include "timeval.h"
   
Line 193 Line 235
 #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 unix  #define AGEGOMP 10. /* Minimal age for Gompertz adjustment */
   #ifdef UNIX
 #define DIRSEPARATOR '/'  #define DIRSEPARATOR '/'
 #define ODIRSEPARATOR '\\'  #define ODIRSEPARATOR '\\'
 #else  #else
Line 204 Line 247
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
   
 char version[]="Imach version 0.96b, June 2003, INED-EUROREVES ";  char version[]="Imach version 0.97c, September 2004, INED-EUROREVES ";
 char fullversion[]="$Revision$ $Date$";   char fullversion[]="$Revision$ $Date$"; 
 int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */  int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
 int nvar;  int nvar;
Line 234  int globpr; /* Global variable for print Line 277  int globpr; /* Global variable for print
 double fretone; /* Only one call to likelihood */  double fretone; /* Only one call to likelihood */
 long ipmx; /* Number of contributions */  long ipmx; /* Number of contributions */
 double sw; /* Sum of weights */  double sw; /* Sum of weights */
   char filerespow[FILENAMELENGTH];
 char fileresilk[FILENAMELENGTH]; /* File of individual contributions to the likelihood */  char fileresilk[FILENAMELENGTH]; /* File of individual contributions to the likelihood */
 FILE *ficresilk;  FILE *ficresilk;
 FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;  FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;
Line 248  char fileresvpl[FILENAMELENGTH]; Line 292  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;
   
Line 296  static double maxarg1,maxarg2; Line 340  static double maxarg1,maxarg2;
 static double sqrarg;  static double sqrarg;
 #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 :sqrarg*sqrarg)  #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 :sqrarg*sqrarg)
 #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}   #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;} 
   int agegomp= AGEGOMP;
   
 int imx;   int imx; 
 int stepm;  int stepm=1;
 /* Stepm, step in month: minimum step interpolation*/  /* Stepm, step in month: minimum step interpolation*/
   
 int estepm;  int estepm;
Line 306  int estepm; Line 351  int estepm;
   
 int m,nb;  int m,nb;
 long *num;  long *num;
 int firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;  int firstpass=0, lastpass=4,*cod, *ncodemax, *Tage,*cens;
 double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;  double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
 double **pmmij, ***probs;  double **pmmij, ***probs;
   double *ageexmed,*agecens;
 double dateintmean=0;  double dateintmean=0;
   
 double *weight;  double *weight;
Line 322  double ftolhess; /* Tolerance for comput Line 368  double ftolhess; /* Tolerance for comput
 /**************** split *************************/  /**************** split *************************/
 static  int split( char *path, char *dirc, char *name, char *ext, char *finame )  static  int split( char *path, char *dirc, char *name, char *ext, char *finame )
 {  {
     /* From a file name with full path (either Unix or Windows) we extract the directory (dirc)
        the name of the file (name), its extension only (ext) and its first part of the name (finame)
     */ 
   char  *ss;                            /* pointer */    char  *ss;                            /* pointer */
   int   l1, l2;                         /* length counters */    int   l1, l2;                         /* length counters */
   
Line 353  static int split( char *path, char *dirc Line 402  static int split( char *path, char *dirc
 #endif  #endif
   */    */
   ss = strrchr( name, '.' );            /* find last / */    ss = strrchr( name, '.' );            /* find last / */
   ss++;    if (ss >0){
   strcpy(ext,ss);                       /* save extension */      ss++;
   l1= strlen( name);      strcpy(ext,ss);                     /* save extension */
   l2= strlen(ss)+1;      l1= strlen( name);
   strncpy( finame, name, l1-l2);      l2= strlen(ss)+1;
   finame[l1-l2]= 0;      strncpy( finame, name, l1-l2);
       finame[l1-l2]= 0;
     }
   return( 0 );                          /* we're done */    return( 0 );                          /* we're done */
 }  }
   
Line 391  int nbocc(char *s, char occ) Line 442  int nbocc(char *s, char occ)
   
 void cutv(char *u,char *v, char*t, char occ)  void cutv(char *u,char *v, char*t, char occ)
 {  {
   /* cuts string t into u and v where u is ended by char occ excluding it    /* cuts string t into u and v where u ends before first occurence of char 'occ' 
      and v is after occ excluding it too : ex cutv(u,v,"abcdef2ghi2j",2)       and v starts after first occurence of char 'occ' : ex cutv(u,v,"abcdef2ghi2j",'2')
      gives u="abcedf" and v="ghi2j" */       gives u="abcedf" and v="ghi2j" */
   int i,lg,j,p=0;    int i,lg,j,p=0;
   i=0;    i=0;
Line 821  void powell(double p[], double **xi, int Line 872  void powell(double p[], double **xi, int
     last_time=curr_time;      last_time=curr_time;
     (void) gettimeofday(&curr_time,&tzp);      (void) gettimeofday(&curr_time,&tzp);
     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);      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(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);      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]);
       fprintf(ficrespow," %.12lf", p[i]);        fprintf(ficrespow," %.12lf", p[i]);
Line 833  void powell(double p[], double **xi, int Line 885  void powell(double p[], double **xi, int
     fprintf(ficrespow,"\n");fflush(ficrespow);      fprintf(ficrespow,"\n");fflush(ficrespow);
     if(*iter <=3){      if(*iter <=3){
       tm = *localtime(&curr_time.tv_sec);        tm = *localtime(&curr_time.tv_sec);
       strcpy(strcurr,asctime(&tmf));        strcpy(strcurr,asctime(&tm));
 /*       asctime_r(&tm,strcurr); */  /*       asctime_r(&tm,strcurr); */
       forecast_time=curr_time;        forecast_time=curr_time; 
       itmp = strlen(strcurr);        itmp = strlen(strcurr);
       if(strcurr[itmp-1]=='\n')        if(strcurr[itmp-1]=='\n')  /* Windows outputs with a new line */
         strcurr[itmp-1]='\0';          strcurr[itmp-1]='\0';
       printf("\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);        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);        fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);
Line 849  void powell(double p[], double **xi, int Line 901  void powell(double p[], double **xi, int
         itmp = strlen(strfor);          itmp = strlen(strfor);
         if(strfor[itmp-1]=='\n')          if(strfor[itmp-1]=='\n')
         strfor[itmp-1]='\0';          strfor[itmp-1]='\0';
         printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s or\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);          printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(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);          fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
       }        }
     }      }
     for (i=1;i<=n;i++) {       for (i=1;i<=n;i++) { 
Line 1015  double **pmij(double **ps, double *cov, Line 1067  double **pmij(double **ps, double *cov,
   int i,j,j1, nc, ii, jj;    int i,j,j1, nc, ii, jj;
   
     for(i=1; i<= nlstate; i++){      for(i=1; i<= nlstate; i++){
     for(j=1; j<i;j++){        for(j=1; j<i;j++){
       for (nc=1, s2=0.;nc <=ncovmodel; nc++){          for (nc=1, s2=0.;nc <=ncovmodel; nc++){
         /*s2 += param[i][j][nc]*cov[nc];*/            /*s2 += param[i][j][nc]*cov[nc];*/
         s2 += x[(i-1)*nlstate*ncovmodel+(j-1)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];            s2 += x[(i-1)*nlstate*ncovmodel+(j-1)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
         /*printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2);*/  /*       printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2); */
       }          }
       ps[i][j]=s2;          ps[i][j]=s2;
       /*printf("s1=%.17e, s2=%.17e\n",s1,s2);*/  /*      printf("s1=%.17e, s2=%.17e\n",s1,s2); */
     }        }
     for(j=i+1; j<=nlstate+ndeath;j++){        for(j=i+1; j<=nlstate+ndeath;j++){
       for (nc=1, s2=0.;nc <=ncovmodel; nc++){          for (nc=1, s2=0.;nc <=ncovmodel; nc++){
         s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];            s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
         /*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/  /*        printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2); */
           }
           ps[i][j]=s2;
       }        }
       ps[i][j]=s2;  
     }      }
   }  
     /*ps[3][2]=1;*/      /*ps[3][2]=1;*/
       
   for(i=1; i<= nlstate; i++){      for(i=1; i<= nlstate; i++){
      s1=0;        s1=0;
     for(j=1; j<i; j++)        for(j=1; j<i; j++)
       s1+=exp(ps[i][j]);          s1+=exp(ps[i][j]);
     for(j=i+1; j<=nlstate+ndeath; j++)        for(j=i+1; j<=nlstate+ndeath; j++)
       s1+=exp(ps[i][j]);          s1+=exp(ps[i][j]);
     ps[i][i]=1./(s1+1.);        ps[i][i]=1./(s1+1.);
     for(j=1; j<i; j++)        for(j=1; j<i; j++)
       ps[i][j]= exp(ps[i][j])*ps[i][i];          ps[i][j]= exp(ps[i][j])*ps[i][i];
     for(j=i+1; j<=nlstate+ndeath; j++)        for(j=i+1; j<=nlstate+ndeath; j++)
       ps[i][j]= exp(ps[i][j])*ps[i][i];          ps[i][j]= exp(ps[i][j])*ps[i][i];
     /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */        /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
   } /* end i */      } /* end i */
       
   for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){      for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
     for(jj=1; jj<= nlstate+ndeath; jj++){        for(jj=1; jj<= nlstate+ndeath; jj++){
       ps[ii][jj]=0;          ps[ii][jj]=0;
       ps[ii][ii]=1;          ps[ii][ii]=1;
         }
     }      }
   }      
   
   
   /*   for(ii=1; ii<= nlstate+ndeath; ii++){  /*        for(ii=1; ii<= nlstate+ndeath; ii++){ */
     for(jj=1; jj<= nlstate+ndeath; jj++){  /*       for(jj=1; jj<= nlstate+ndeath; jj++){ */
      printf("%lf ",ps[ii][jj]);  /*         printf("ddd %lf ",ps[ii][jj]); */
    }  /*       } */
     printf("\n ");  /*       printf("\n "); */
     }  /*        } */
     printf("\n ");printf("%lf ",cov[2]);*/  /*        printf("\n ");printf("%lf ",cov[2]); */
 /*         /*
   for(i=1; i<= npar; i++) printf("%f ",x[i]);        for(i=1; i<= npar; i++) printf("%f ",x[i]);
   goto end;*/        goto end;*/
     return ps;      return ps;
 }  }
   
Line 1189  double func( double *x) Line 1241  double func( double *x)
         } /* end mult */          } /* end mult */
               
         /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */          /*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.          /* But now since version 0.9 we anticipate for bias at large stepm.
          * If stepm is larger than one month (smallest stepm) and if the exact delay            * 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            * (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           * 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           * 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           * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the
          * probability in order to take into account the bias as a fraction of the way           * 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           * from savm to out if bh is negative or even beyond if bh is positive. bh varies
          * -stepm/2 to stepm/2 .           * -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 the same as for previous versions of Imach.
          * For stepm > 1 the results are less biased than in previous versions.            * For stepm > 1 the results are less biased than in previous versions. 
Line 1204  double func( double *x) Line 1256  double func( double *x)
         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          /* bias bh is positive if real duration
          * is higher than the multiple of stepm and negative otherwise.           * is higher than the multiple of stepm and negative otherwise.
          */           */
         /* 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.+bbh)*out[s1][s2]));*/
         if( s2 > nlstate){           if( s2 > nlstate){ 
           /* i.e. if s2 is a death state and if the date of death is known then the contribution            /* 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                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                step unit time, which is also equal to probability to die before dh 
              and probability to die before dh-stepm .                minus probability to die before dh-stepm . 
              In version up to 0.92 likelihood was computed               In version up to 0.92 likelihood was computed
         as if date of death was unknown. Death was treated as any other          as if date of death was unknown. Death was treated as any other
         health state: the date of the interview describes the actual state          health state: the date of the interview describes the actual state
Line 1508  void mlikeli(FILE *ficres,double p[], in Line 1560  void mlikeli(FILE *ficres,double p[], in
   double **xi;    double **xi;
   double fret;    double fret;
   double fretone; /* Only one call to likelihood */    double fretone; /* Only one call to likelihood */
   char filerespow[FILENAMELENGTH];    /*  char filerespow[FILENAMELENGTH];*/
   xi=matrix(1,npar,1,npar);    xi=matrix(1,npar,1,npar);
   for (i=1;i<=npar;i++)    for (i=1;i<=npar;i++)
     for (j=1;j<=npar;j++)      for (j=1;j<=npar;j++)
Line 1543  void hesscov(double **matcov, double p[] Line 1595  void hesscov(double **matcov, double p[]
   int i, j,jk;    int i, j,jk;
   int *indx;    int *indx;
   
   double hessii(double p[], double delta, int theta, double delti[]);    double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar);
   double hessij(double p[], double delti[], int i, int j);    double hessij(double p[], double delti[], int i, int j,double (*func)(double []),int npar);
   void lubksb(double **a, int npar, int *indx, double b[]) ;    void lubksb(double **a, int npar, int *indx, double b[]) ;
   void ludcmp(double **a, int npar, int *indx, double *d) ;    void ludcmp(double **a, int npar, int *indx, double *d) ;
     double gompertz(double p[]);
   hess=matrix(1,npar,1,npar);    hess=matrix(1,npar,1,npar);
   
   printf("\nCalculation of the hessian matrix. Wait...\n");    printf("\nCalculation of the hessian matrix. Wait...\n");
Line 1555  void hesscov(double **matcov, double p[] Line 1607  void hesscov(double **matcov, double p[]
   for (i=1;i<=npar;i++){    for (i=1;i<=npar;i++){
     printf("%d",i);fflush(stdout);      printf("%d",i);fflush(stdout);
     fprintf(ficlog,"%d",i);fflush(ficlog);      fprintf(ficlog,"%d",i);fflush(ficlog);
     hess[i][i]=hessii(p,ftolhess,i,delti);     
     /*printf(" %f ",p[i]);*/       hess[i][i]=hessii(p,ftolhess,i,delti,func,npar);
     /*printf(" %lf ",hess[i][i]);*/      
       /*  printf(" %f ",p[i]);
           printf(" %lf %lf %lf",hess[i][i],ftolhess,delti[i]);*/
   }    }
       
   for (i=1;i<=npar;i++) {    for (i=1;i<=npar;i++) {
Line 1565  void hesscov(double **matcov, double p[] Line 1619  void hesscov(double **matcov, double p[]
       if (j>i) {         if (j>i) { 
         printf(".%d%d",i,j);fflush(stdout);          printf(".%d%d",i,j);fflush(stdout);
         fprintf(ficlog,".%d%d",i,j);fflush(ficlog);          fprintf(ficlog,".%d%d",i,j);fflush(ficlog);
         hess[i][j]=hessij(p,delti,i,j);          hess[i][j]=hessij(p,delti,i,j,func,npar);
           
         hess[j][i]=hess[i][j];              hess[j][i]=hess[i][j];    
         /*printf(" %lf ",hess[i][j]);*/          /*printf(" %lf ",hess[i][j]);*/
       }        }
Line 1636  void hesscov(double **matcov, double p[] Line 1691  void hesscov(double **matcov, double p[]
 }  }
   
 /*************** hessian matrix ****************/  /*************** hessian matrix ****************/
 double hessii( double x[], double delta, int theta, double delti[])  double hessii(double x[], double delta, int theta, double delti[], double (*func)(double []), int npar)
 {  {
   int i;    int i;
   int l=1, lmax=20;    int l=1, lmax=20;
   double k1,k2;    double k1,k2;
   double p2[NPARMAX+1];    double p2[NPARMAX+1];
   double res;    double res;
   double delt, delts, nkhi=10.,nkhif=1., khi=1.e-4;    double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4;
   double fx;    double fx;
   int k=0,kmax=10;    int k=0,kmax=10;
   double l1;    double l1;
Line 1683  double hessii( double x[], double delta, Line 1738  double hessii( double x[], double delta,
       
 }  }
   
 double hessij( double x[], double delti[], int thetai,int thetaj)  double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar)
 {  {
   int i;    int i;
   int l=1, l1, lmax=20;    int l=1, l1, lmax=20;
Line 2887  To be simple, these graphs help to under Line 2942  To be simple, these graphs help to under
         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 3098  void printinghtml(char fileres[], char t Line 3153  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 3146  fprintf(fichtm," \n<ul><li><b>Graphs</b> Line 3199  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 3176  health expectancies in states (1) and (2 Line 3233  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 3196  fprintf(fichtm," <ul><li><b>Graphs</b></ Line 3254  fprintf(fichtm," <ul><li><b>Graphs</b></
 prevalence (with 95%% confidence 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,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),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 3853  void prwizard(int ncovmodel, int nlstate Line 3914  void prwizard(int ncovmodel, int nlstate
   } /* end itimes */    } /* end itimes */
   
 } /* end of prwizard */  } /* end of prwizard */
   /******************* Gompertz Likelihood ******************************/
   double gompertz(double x[])
   { 
     double A,B,L=0.0,sump=0.,num=0.;
     int i,n=0; /* n is the size of the sample */
     for (i=0;i<=imx-1 ; i++) {
       sump=sump+weight[i];
       sump=sump+1;
       num=num+1;
     }
    
    
     /* for (i=1; i<=imx; i++) 
        if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/
   
     for (i=0;i<=imx-1 ; i++)
       {
         if (cens[i]==1 & wav[i]>1)
           A=-x[1]/(x[2])*
             (exp(x[2]/YEARM*(agecens[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12)));
         
         if (cens[i]==0 & wav[i]>1)
           A=-x[1]/(x[2])*
                (exp(x[2]/YEARM*(agedc[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12)))
             +log(x[1]/YEARM)+x[2]/YEARM*(agedc[i]*12-agegomp*12)+log(YEARM);      
         
         if (wav[i]>1 & agecens[i]>15) {
           L=L+A*weight[i];
           /*      printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/
         }
       }
   
    /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/
    
     return -2*L*num/sump;
   }
   
   /******************* Printing html file ***********/
   void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \
                     int lastpass, int stepm, int weightopt, char model[],\
                     int imx,  double p[],double **matcov){
     int i;
   
     fprintf(fichtm,"<ul><li><h4>Result files </h4>\n Force of mortality. Parameters of the Gompertz fit (with confidence interval in brackets):<br>");
     fprintf(fichtm,"  mu(age) =%lf*exp(%lf*(age-%d)) per year<br><br>",p[1],p[2],agegomp);
     for (i=1;i<=2;i++) 
       fprintf(fichtm," p[%d] = %lf [%f ; %f]<br>\n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
     fprintf(fichtm,"<br><br><img src=\"graphmort.png\">");
     fprintf(fichtm,"</ul>");
     fflush(fichtm);
   }
   
   /******************* Gnuplot file **************/
   void printinggnuplotmort(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
   
     char dirfileres[132],optfileres[132];
     int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
     int ng;
   
   
     /*#ifdef windows */
     fprintf(ficgp,"cd \"%s\" \n",pathc);
       /*#endif */
   
   
     strcpy(dirfileres,optionfilefiname);
     strcpy(optfileres,"vpl");
     fprintf(ficgp,"set out \"graphmort.png\"\n "); 
     fprintf(ficgp,"set xlabel \"Age\"\n set ylabel \"Force of mortality (per year)\" \n "); 
     fprintf(ficgp, "set ter png small\n set log y\n"); 
     fprintf(ficgp, "set size 0.65,0.65\n");
     fprintf(ficgp,"plot [%d:100] %lf*exp(%lf*(x-%d))",agegomp,p[1],p[2],agegomp);
   
   } 
   
   
   
   
 /***********************************************/  /***********************************************/
Line 3866  int main(int argc, char *argv[]) Line 4003  int main(int argc, char *argv[])
   int jj, ll, li, lj, lk, 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;    int itimes;
     int NDIM=2;
   
   char ca[32], cb[32], cc[32];    char ca[32], cb[32], cc[32];
   /*  FILE *fichtm; *//* Html File */    /*  FILE *fichtm; *//* Html File */
Line 3882  int main(int argc, char *argv[]) Line 4020  int main(int argc, char *argv[])
   int *indx;    int *indx;
   char line[MAXLINE], linepar[MAXLINE];    char line[MAXLINE], linepar[MAXLINE];
   char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];    char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];
   char pathr[MAXLINE];     char pathr[MAXLINE], pathimach[MAXLINE]; 
   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 */
   int c,  h , cpt,l;    int c,  h , cpt,l;
Line 3909  int main(int argc, char *argv[]) Line 4047  int main(int argc, char *argv[])
   double *epj, vepp;    double *epj, vepp;
   double kk1, kk2;    double kk1, kk2;
   double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;    double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
     double **ximort;
   char *alph[]={"a","a","b","c","d","e"}, str[4];    char *alph[]={"a","a","b","c","d","e"}, str[4];
     int *dcwave;
   
   char z[1]="c", occ;    char z[1]="c", occ;
   
Line 3968  int main(int argc, char *argv[]) Line 4106  int main(int argc, char *argv[])
     printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/      printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
   /* cutv(path,optionfile,pathtot,'\\');*/    /* cutv(path,optionfile,pathtot,'\\');*/
   
     split(argv[0],pathimach,optionfile,optionfilext,optionfilefiname);
    /*   strcpy(pathimach,argv[0]); */
   split(pathtot,path,optionfile,optionfilext,optionfilefiname);    split(pathtot,path,optionfile,optionfilext,optionfilefiname);
   printf("pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);    printf("pathimach=%s, pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
   chdir(path);    chdir(path);
   strcpy(command,"mkdir ");    strcpy(command,"mkdir ");
   strcat(command,optionfilefiname);    strcat(command,optionfilefiname);
Line 3994  int main(int argc, char *argv[]) Line 4134  int main(int argc, char *argv[])
   }    }
   fprintf(ficlog,"Log filename:%s\n",filelog);    fprintf(ficlog,"Log filename:%s\n",filelog);
   fprintf(ficlog,"\n%s\n%s",version,fullversion);    fprintf(ficlog,"\n%s\n%s",version,fullversion);
   fprintf(ficlog,"\nEnter the parameter file name: ");    fprintf(ficlog,"\nEnter the parameter file name: \n");
   fprintf(ficlog,"pathtot=%s\n\    fprintf(ficlog,"pathimach=%s\npathtot=%s\n\
  path=%s \n\   path=%s \n\
  optionfile=%s\n\   optionfile=%s\n\
  optionfilext=%s\n\   optionfilext=%s\n\
  optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);   optionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
   
   printf("Local time (at start):%s",strstart);    printf("Local time (at start):%s",strstart);
   fprintf(ficlog,"Local time (at start): %s",strstart);    fprintf(ficlog,"Local time (at start): %s",strstart);
Line 4067  int main(int argc, char *argv[]) Line 4207  int main(int argc, char *argv[])
   
   ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */    ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */
   nvar=ncovmodel-1; /* Suppressing age as a basic covariate */    nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
      npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/
   
     delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
     delti=delti3[1][1];
     /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/
   if(mle==-1){ /* Print a wizard for help writing covariance matrix */    if(mle==-1){ /* Print a wizard for help writing covariance matrix */
     prwizard(ncovmodel, nlstate, ndeath, model, ficparo);      prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
     printf(" You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso);      printf(" You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
     fprintf(ficlog," You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso);      fprintf(ficlog," You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
       free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); 
     fclose (ficparo);      fclose (ficparo);
     fclose (ficlog);      fclose (ficlog);
     exit(0);      exit(0);
   }    }
   /* Read guess parameters */    else if(mle==-3) {
   /* Reads comments: lines beginning with '#' */      prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
   while((c=getc(ficpar))=='#' && c!= EOF){      printf(" You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
     ungetc(c,ficpar);      fprintf(ficlog," You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
     fgets(line, MAXLINE, ficpar);      param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
     numlinepar++;      matcov=matrix(1,npar,1,npar);
     puts(line);  
     fputs(line,ficparo);  
     fputs(line,ficlog);  
   }    }
   ungetc(c,ficpar);    else{
       /* Read guess parameters */
   param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);      /* Reads comments: lines beginning with '#' */
   for(i=1; i <=nlstate; i++){      while((c=getc(ficpar))=='#' && c!= EOF){
     j=0;        ungetc(c,ficpar);
     for(jj=1; jj <=nlstate+ndeath; jj++){        fgets(line, MAXLINE, ficpar);
       if(jj==i) continue;  
       j++;  
       fscanf(ficpar,"%1d%1d",&i1,&j1);  
       if ((i1 != i) && (j1 != j)){  
         printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);  
         exit(1);  
       }  
       fprintf(ficparo,"%1d%1d",i1,j1);  
       if(mle==1)  
         printf("%1d%1d",i,j);  
       fprintf(ficlog,"%1d%1d",i,j);  
       for(k=1; k<=ncovmodel;k++){  
         fscanf(ficpar," %lf",&param[i][j][k]);  
         if(mle==1){  
           printf(" %lf",param[i][j][k]);  
           fprintf(ficlog," %lf",param[i][j][k]);  
         }  
         else  
           fprintf(ficlog," %lf",param[i][j][k]);  
         fprintf(ficparo," %lf",param[i][j][k]);  
       }  
       fscanf(ficpar,"\n");  
       numlinepar++;        numlinepar++;
       if(mle==1)        puts(line);
         printf("\n");        fputs(line,ficparo);
       fprintf(ficlog,"\n");        fputs(line,ficlog);
       fprintf(ficparo,"\n");  
     }      }
   }        ungetc(c,ficpar);
   fflush(ficlog);      
       param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
       for(i=1; i <=nlstate; i++){
         j=0;
         for(jj=1; jj <=nlstate+ndeath; jj++){
           if(jj==i) continue;
           j++;
           fscanf(ficpar,"%1d%1d",&i1,&j1);
           if ((i1 != i) && (j1 != j)){
             printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);
             exit(1);
           }
           fprintf(ficparo,"%1d%1d",i1,j1);
           if(mle==1)
             printf("%1d%1d",i,j);
           fprintf(ficlog,"%1d%1d",i,j);
           for(k=1; k<=ncovmodel;k++){
             fscanf(ficpar," %lf",&param[i][j][k]);
             if(mle==1){
               printf(" %lf",param[i][j][k]);
               fprintf(ficlog," %lf",param[i][j][k]);
             }
             else
               fprintf(ficlog," %lf",param[i][j][k]);
             fprintf(ficparo," %lf",param[i][j][k]);
           }
           fscanf(ficpar,"\n");
           numlinepar++;
           if(mle==1)
             printf("\n");
           fprintf(ficlog,"\n");
           fprintf(ficparo,"\n");
         }
       }  
       fflush(ficlog);
   
   npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/  
   
   p=param[1][1];      p=param[1][1];
         
   /* Reads comments: lines beginning with '#' */      /* Reads comments: lines beginning with '#' */
   while((c=getc(ficpar))=='#' && c!= EOF){      while((c=getc(ficpar))=='#' && c!= EOF){
         ungetc(c,ficpar);
         fgets(line, MAXLINE, ficpar);
         numlinepar++;
         puts(line);
         fputs(line,ficparo);
         fputs(line,ficlog);
       }
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);  
     numlinepar++;  
     puts(line);  
     fputs(line,ficparo);  
     fputs(line,ficlog);  
   }  
   ungetc(c,ficpar);  
   
   delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);      for(i=1; i <=nlstate; i++){
   /* delti=vector(1,npar); *//* Scale of each paramater (output from hesscov) */        for(j=1; j <=nlstate+ndeath-1; j++){
   for(i=1; i <=nlstate; i++){          fscanf(ficpar,"%1d%1d",&i1,&j1);
     for(j=1; j <=nlstate+ndeath-1; j++){          if ((i1-i)*(j1-j)!=0){
       fscanf(ficpar,"%1d%1d",&i1,&j1);            printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);
       if ((i1-i)*(j1-j)!=0){            exit(1);
         printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);          }
         exit(1);          printf("%1d%1d",i,j);
       }          fprintf(ficparo,"%1d%1d",i1,j1);
       printf("%1d%1d",i,j);          fprintf(ficlog,"%1d%1d",i1,j1);
       fprintf(ficparo,"%1d%1d",i1,j1);          for(k=1; k<=ncovmodel;k++){
       fprintf(ficlog,"%1d%1d",i1,j1);            fscanf(ficpar,"%le",&delti3[i][j][k]);
       for(k=1; k<=ncovmodel;k++){            printf(" %le",delti3[i][j][k]);
         fscanf(ficpar,"%le",&delti3[i][j][k]);            fprintf(ficparo," %le",delti3[i][j][k]);
         printf(" %le",delti3[i][j][k]);            fprintf(ficlog," %le",delti3[i][j][k]);
         fprintf(ficparo," %le",delti3[i][j][k]);          }
         fprintf(ficlog," %le",delti3[i][j][k]);          fscanf(ficpar,"\n");
           numlinepar++;
           printf("\n");
           fprintf(ficparo,"\n");
           fprintf(ficlog,"\n");
       }        }
       fscanf(ficpar,"\n");  
       numlinepar++;  
       printf("\n");  
       fprintf(ficparo,"\n");  
       fprintf(ficlog,"\n");  
     }      }
   }      fflush(ficlog);
   fflush(ficlog);  
   
   delti=delti3[1][1];      delti=delti3[1][1];
   
   
   /* free_ma3x(delti3,1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */ /* Hasn't to to freed here otherwise delti is no more allocated */      /* free_ma3x(delti3,1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */ /* Hasn't to to freed here otherwise delti is no more allocated */
       
   /* Reads comments: lines beginning with '#' */      /* Reads comments: lines beginning with '#' */
   while((c=getc(ficpar))=='#' && c!= EOF){      while((c=getc(ficpar))=='#' && c!= EOF){
         ungetc(c,ficpar);
         fgets(line, MAXLINE, ficpar);
         numlinepar++;
         puts(line);
         fputs(line,ficparo);
         fputs(line,ficlog);
       }
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);  
     numlinepar++;  
     puts(line);  
     fputs(line,ficparo);  
     fputs(line,ficlog);  
   }  
   ungetc(c,ficpar);  
       
   matcov=matrix(1,npar,1,npar);      matcov=matrix(1,npar,1,npar);
   for(i=1; i <=npar; i++){      for(i=1; i <=npar; i++){
     fscanf(ficpar,"%s",&str);        fscanf(ficpar,"%s",&str);
     if(mle==1)        if(mle==1)
       printf("%s",str);          printf("%s",str);
     fprintf(ficlog,"%s",str);        fprintf(ficlog,"%s",str);
     fprintf(ficparo,"%s",str);        fprintf(ficparo,"%s",str);
     for(j=1; j <=i; j++){        for(j=1; j <=i; j++){
       fscanf(ficpar," %le",&matcov[i][j]);          fscanf(ficpar," %le",&matcov[i][j]);
       if(mle==1){          if(mle==1){
         printf(" %.5le",matcov[i][j]);            printf(" %.5le",matcov[i][j]);
           }
           fprintf(ficlog," %.5le",matcov[i][j]);
           fprintf(ficparo," %.5le",matcov[i][j]);
       }        }
       fprintf(ficlog," %.5le",matcov[i][j]);        fscanf(ficpar,"\n");
       fprintf(ficparo," %.5le",matcov[i][j]);        numlinepar++;
         if(mle==1)
           printf("\n");
         fprintf(ficlog,"\n");
         fprintf(ficparo,"\n");
     }      }
     fscanf(ficpar,"\n");      for(i=1; i <=npar; i++)
     numlinepar++;        for(j=i+1;j<=npar;j++)
           matcov[i][j]=matcov[j][i];
       
     if(mle==1)      if(mle==1)
       printf("\n");        printf("\n");
     fprintf(ficlog,"\n");      fprintf(ficlog,"\n");
     fprintf(ficparo,"\n");  
   }  
   for(i=1; i <=npar; i++)  
     for(j=i+1;j<=npar;j++)  
       matcov[i][j]=matcov[j][i];  
      
   if(mle==1)  
     printf("\n");  
   fprintf(ficlog,"\n");  
   
   fflush(ficlog);  
   
   /*-------- Rewriting paramater file ----------*/  
   strcpy(rfileres,"r");    /* "Rparameterfile */  
   strcat(rfileres,optionfilefiname);    /* Parameter file first name*/  
   strcat(rfileres,".");    /* */  
   strcat(rfileres,optionfilext);    /* Other files have txt extension */  
   if((ficres =fopen(rfileres,"w"))==NULL) {  
     printf("Problem writing new parameter file: %s\n", fileres);goto end;  
     fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;  
   }  
   fprintf(ficres,"#%s\n",version);  
           
       fflush(ficlog);
       
       /*-------- Rewriting parameter file ----------*/
       strcpy(rfileres,"r");    /* "Rparameterfile */
       strcat(rfileres,optionfilefiname);    /* Parameter file first name*/
       strcat(rfileres,".");    /* */
       strcat(rfileres,optionfilext);    /* Other files have txt extension */
       if((ficres =fopen(rfileres,"w"))==NULL) {
         printf("Problem writing new parameter file: %s\n", fileres);goto end;
         fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;
       }
       fprintf(ficres,"#%s\n",version);
     }    /* End of mle != -3 */
   
   /*-------- data file ----------*/    /*-------- data file ----------*/
   if((fic=fopen(datafile,"r"))==NULL)    {    if((fic=fopen(datafile,"r"))==NULL)    {
     printf("Problem with datafile: %s\n", datafile);goto end;      printf("Problem with datafile: %s\n", datafile);goto end;
Line 4251  int main(int argc, char *argv[]) Line 4402  int main(int argc, char *argv[])
   i=1;    i=1;
   while (fgets(line, MAXLINE, fic) != NULL)    {    while (fgets(line, MAXLINE, fic) != NULL)    {
     if ((i >= firstobs) && (i <=lastobs)) {      if ((i >= firstobs) && (i <=lastobs)) {
                 for(j=0; line[j] != '\n';j++){  /* Untabifies line */
           if(line[j] == '\t')
             line[j] = ' ';
         }
       for (j=maxwav;j>=1;j--){        for (j=maxwav;j>=1;j--){
         cutv(stra, strb,line,' '); s[j][i]=atoi(strb);           cutv(stra, strb,line,' '); s[j][i]=atoi(strb); 
         strcpy(line,stra);          strcpy(line,stra);
Line 4474  int main(int argc, char *argv[]) Line 4628  int main(int argc, char *argv[])
   
 }*/  }*/
   
   
   printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);    printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
   fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);     fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); 
   
     agegomp=(int)agemin;
   free_vector(severity,1,maxwav);    free_vector(severity,1,maxwav);
   free_imatrix(outcome,1,maxwav+1,1,n);    free_imatrix(outcome,1,maxwav+1,1,n);
   free_vector(moisnais,1,n);    free_vector(moisnais,1,n);
Line 4530  int main(int argc, char *argv[]) Line 4686  int main(int argc, char *argv[])
           
   /*------------ gnuplot -------------*/    /*------------ gnuplot -------------*/
   strcpy(optionfilegnuplot,optionfilefiname);    strcpy(optionfilegnuplot,optionfilefiname);
     if(mle==-3)
       strcat(optionfilegnuplot,"-mort");
   strcat(optionfilegnuplot,".gp");    strcat(optionfilegnuplot,".gp");
   
   if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {    if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
     printf("Problem with file %s",optionfilegnuplot);      printf("Problem with file %s",optionfilegnuplot);
   }    }
Line 4543  int main(int argc, char *argv[]) Line 4702  int main(int argc, char *argv[])
   /*--------- index.htm --------*/    /*--------- index.htm --------*/
   
   strcpy(optionfilehtm,optionfilefiname); /* Main html file */    strcpy(optionfilehtm,optionfilefiname); /* Main html file */
     if(mle==-3)
       strcat(optionfilehtm,"-mort");
   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);
Line 4600  Interval (in months) between two waves: Line 4761  Interval (in months) between two waves:
   p=param[1][1]; /* *(*(*(param +1)+1)+0) */    p=param[1][1]; /* *(*(*(param +1)+1)+0) */
   
   globpr=0; /* To get the number ipmx of contributions and the sum of weights*/    globpr=0; /* To get the number ipmx of contributions and the sum of weights*/
   likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */    if (mle==-3){
   printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);      ximort=matrix(1,NDIM,1,NDIM);
   for (k=1; k<=npar;k++)      cens=ivector(1,n);
     printf(" %d %8.5f",k,p[k]);      ageexmed=vector(1,n);
   printf("\n");      agecens=vector(1,n);
   globpr=1; /* to print the contributions */      dcwave=ivector(1,n);
   likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */   
   printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);      for (i=1; i<=imx; i++){
   for (k=1; k<=npar;k++)        dcwave[i]=-1;
     printf(" %d %8.5f",k,p[k]);        for (j=1; j<=lastpass; j++)
   printf("\n");          if (s[j][i]>nlstate) {
   if(mle>=1){ /* Could be 1 or 2 */            dcwave[i]=j;
     mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);            /*    printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/
             break;
           }
       }
   
       for (i=1; i<=imx; i++) {
         if (wav[i]>0){
           ageexmed[i]=agev[mw[1][i]][i];
           j=wav[i];agecens[i]=1.; 
           if (ageexmed[i]>1 & wav[i]>0) agecens[i]=agev[mw[j][i]][i];
           cens[i]=1;
           
           if (ageexmed[i]<1) cens[i]=-1;
           if (agedc[i]< AGESUP & agedc[i]>1 & dcwave[i]>firstpass & dcwave[i]<=lastpass) cens[i]=0 ;
         }
         else cens[i]=-1;
       }
       
       for (i=1;i<=NDIM;i++) {
         for (j=1;j<=NDIM;j++)
           ximort[i][j]=(i == j ? 1.0 : 0.0);
       }
   
       p[1]=0.1; p[2]=0.1;
       /*printf("%lf %lf", p[1], p[2]);*/
       
       
     printf("Powell\n");  fprintf(ficlog,"Powell\n");
     strcpy(filerespow,"pow-mort"); 
     strcat(filerespow,fileres);
     if((ficrespow=fopen(filerespow,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", filerespow);
       fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
   }    }
     fprintf(ficrespow,"# Powell\n# iter -2*LL");
     /*  for (i=1;i<=nlstate;i++)
       for(j=1;j<=nlstate+ndeath;j++)
         if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
     */
     fprintf(ficrespow,"\n");
   
       powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);
       fclose(ficrespow);
           
   /*--------- results files --------------*/      hesscov(matcov, p, NDIM,delti, 1e-4, gompertz); 
   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,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");      for(i=1; i <=NDIM; i++)
   printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");        for(j=i+1;j<=NDIM;j++)
   fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");          matcov[i][j]=matcov[j][i];
   for(i=1,jk=1; i <=nlstate; i++){      
     for(k=1; k <=(nlstate+ndeath); k++){      printf("\nCovariance matrix\n ");
       if (k != i) {      for(i=1; i <=NDIM; i++) {
         printf("%d%d ",i,k);        for(j=1;j<=NDIM;j++){ 
         fprintf(ficlog,"%d%d ",i,k);          printf("%f ",matcov[i][j]);
         fprintf(ficres,"%1d%1d ",i,k);  
         for(j=1; j <=ncovmodel; j++){  
           printf("%f ",p[jk]);  
           fprintf(ficlog,"%f ",p[jk]);  
           fprintf(ficres,"%f ",p[jk]);  
           jk++;   
         }  
         printf("\n");  
         fprintf(ficlog,"\n");  
         fprintf(ficres,"\n");  
       }        }
         printf("\n ");
     }      }
   }      
   if(mle!=0){      printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp);
     /* Computing hessian and covariance matrix */      for (i=1;i<=NDIM;i++) 
     ftolhess=ftol; /* Usually correct */        printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
     hesscov(matcov, p, npar, delti, ftolhess, func);      replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */
   }      printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
   fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");      
   printf("# Scales (for hessian or gradient estimation)\n");      printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \
   fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");                       stepm, weightopt,\
   for(i=1,jk=1; i <=nlstate; i++){                       model,imx,p,matcov);
     for(j=1; j <=nlstate+ndeath; j++){    } /* Endof if mle==-3 */
       if (j!=i) {  
         fprintf(ficres,"%1d%1d",i,j);    else{ /* For mle >=1 */
         printf("%1d%1d",i,j);    
         fprintf(ficlog,"%1d%1d",i,j);      likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
         for(k=1; k<=ncovmodel;k++){      printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
           printf(" %.5e",delti[jk]);      for (k=1; k<=npar;k++)
           fprintf(ficlog," %.5e",delti[jk]);        printf(" %d %8.5f",k,p[k]);
           fprintf(ficres," %.5e",delti[jk]);      printf("\n");
           jk++;      globpr=1; /* to print the contributions */
       likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
       printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
       for (k=1; k<=npar;k++)
         printf(" %d %8.5f",k,p[k]);
       printf("\n");
       if(mle>=1){ /* Could be 1 or 2 */
         mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
       }
       
       /*--------- results files --------------*/
       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,"# 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");
       for(i=1,jk=1; i <=nlstate; i++){
         for(k=1; k <=(nlstate+ndeath); k++){
           if (k != i) {
             printf("%d%d ",i,k);
             fprintf(ficlog,"%d%d ",i,k);
             fprintf(ficres,"%1d%1d ",i,k);
             for(j=1; j <=ncovmodel; j++){
               printf("%f ",p[jk]);
               fprintf(ficlog,"%f ",p[jk]);
               fprintf(ficres,"%f ",p[jk]);
               jk++; 
             }
             printf("\n");
             fprintf(ficlog,"\n");
             fprintf(ficres,"\n");
         }          }
         printf("\n");  
         fprintf(ficlog,"\n");  
         fprintf(ficres,"\n");  
       }        }
     }      }
   }      if(mle!=0){
            /* Computing hessian and covariance matrix */
   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");        ftolhess=ftol; /* Usually correct */
   if(mle>=1)        hesscov(matcov, p, npar, delti, ftolhess, func);
     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(ficres,"# Scales (for hessian or gradient estimation)\n");
 /* # 121 Var(a12)\n\ */      printf("# Scales (for hessian or gradient estimation)\n");
 /* # 122 Cov(b12,a12) Var(b12)\n\ */      fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");
 /* # 131 Cov(a13,a12) Cov(a13,b12, Var(a13)\n\ */      for(i=1,jk=1; i <=nlstate; i++){
 /* # 132 Cov(b13,a12) Cov(b13,b12, Cov(b13,a13) Var(b13)\n\ */  
 /* # 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\ */  
 /* # 232 Cov(a23,a12) Cov(a23,b12, Cov(a23,a13) Cov(a23,b13) Cov(a23,a21) Cov(a23,b21) Var(a23)\n\ */  
 /* # 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n" */  
   
   
 /* Just to have a covariance matrix which will be more understandable  
    even is we still don't want to manage dictionary of variables  
 */  
   for(itimes=1;itimes<=2;itimes++){  
     jj=0;  
     for(i=1; i <=nlstate; i++){  
       for(j=1; j <=nlstate+ndeath; j++){        for(j=1; j <=nlstate+ndeath; j++){
         if(j==i) continue;          if (j!=i) {
         for(k=1; k<=ncovmodel;k++){            fprintf(ficres,"%1d%1d",i,j);
           jj++;            printf("%1d%1d",i,j);
           ca[0]= k+'a'-1;ca[1]='\0';            fprintf(ficlog,"%1d%1d",i,j);
           if(itimes==1){            for(k=1; k<=ncovmodel;k++){
             if(mle>=1)              printf(" %.5e",delti[jk]);
               printf("#%1d%1d%d",i,j,k);              fprintf(ficlog," %.5e",delti[jk]);
             fprintf(ficlog,"#%1d%1d%d",i,j,k);              fprintf(ficres," %.5e",delti[jk]);
             fprintf(ficres,"#%1d%1d%d",i,j,k);              jk++;
           }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;            printf("\n");
           for(li=1;li <=nlstate; li++){            fprintf(ficlog,"\n");
             for(lj=1;lj <=nlstate+ndeath; lj++){            fprintf(ficres,"\n");
               if(lj==li) continue;          }
               for(lk=1;lk<=ncovmodel;lk++){        }
                 ll++;      }
                 if(ll<=jj){      
                   cb[0]= lk +'a'-1;cb[1]='\0';      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(ll<jj){      if(mle>=1)
                     if(itimes==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");
                       if(mle>=1)      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");
                         printf(" Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj);      /* # 121 Var(a12)\n\ */
                       fprintf(ficlog," Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj);      /* # 122 Cov(b12,a12) Var(b12)\n\ */
                       fprintf(ficres," Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj);      /* # 131 Cov(a13,a12) Cov(a13,b12, Var(a13)\n\ */
                     }else{      /* # 132 Cov(b13,a12) Cov(b13,b12, Cov(b13,a13) Var(b13)\n\ */
                       if(mle>=1)      /* # 212 Cov(a21,a12) Cov(a21,b12, Cov(a21,a13) Cov(a21,b13) Var(a21)\n\ */
                         printf(" %.5e",matcov[jj][ll]);       /* # 212 Cov(b21,a12) Cov(b21,b12, Cov(b21,a13) Cov(b21,b13) Cov(b21,a21) Var(b21)\n\ */
                       fprintf(ficlog," %.5e",matcov[jj][ll]);       /* # 232 Cov(a23,a12) Cov(a23,b12, Cov(a23,a13) Cov(a23,b13) Cov(a23,a21) Cov(a23,b21) Var(a23)\n\ */
                       fprintf(ficres," %.5e",matcov[jj][ll]);       /* # 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n" */
                     }      
                   }else{      
                     if(itimes==1){      /* Just to have a covariance matrix which will be more understandable
                       if(mle>=1)         even is we still don't want to manage dictionary of variables
                         printf(" Var(%s%1d%1d)",ca,i,j);      */
                       fprintf(ficlog," Var(%s%1d%1d)",ca,i,j);      for(itimes=1;itimes<=2;itimes++){
                       fprintf(ficres," Var(%s%1d%1d)",ca,i,j);        jj=0;
         for(i=1; i <=nlstate; i++){
           for(j=1; j <=nlstate+ndeath; j++){
             if(j==i) continue;
             for(k=1; k<=ncovmodel;k++){
               jj++;
               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{                      }else{
                       if(mle>=1)                        if(itimes==1){
                         printf(" %.5e",matcov[jj][ll]);                           if(mle>=1)
                       fprintf(ficlog," %.5e",matcov[jj][ll]);                             printf(" Var(%s%1d%1d)",ca,i,j);
                       fprintf(ficres," %.5e",matcov[jj][ll]);                           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 lk */                } /* end lj */
             } /* end lj */              } /* end li */
           } /* end li */              if(mle>=1)
           if(mle>=1)                printf("\n");
             printf("\n");              fprintf(ficlog,"\n");
           fprintf(ficlog,"\n");              fprintf(ficres,"\n");
           fprintf(ficres,"\n");              numlinepar++;
           numlinepar++;            } /* end k*/
         } /* end k*/          } /*end j */
       } /*end j */        } /* end i */
     } /* end i */      } /* end itimes */
   } /* end itimes */      
       fflush(ficlog);
   fflush(ficlog);      fflush(ficres);
   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);        puts(line);
     puts(line);        fputs(line,ficparo);
     fputs(line,ficparo);      }
   }  
   ungetc(c,ficpar);  
   
   estepm=0;  
   fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);  
   if (estepm==0 || estepm < stepm) estepm=stepm;  
   if (fage <= 2) {  
     bage = ageminpar;  
     fage = agemaxpar;  
   }  
      
   fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");  
   fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);  
   fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);  
      
   while((c=getc(ficpar))=='#' && c!= EOF){  
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      
     puts(line);      estepm=0;
     fputs(line,ficparo);      fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
   }      if (estepm==0 || estepm < stepm) estepm=stepm;
   ungetc(c,ficpar);      if (fage <= 2) {
           bage = ageminpar;
   fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);        fage = agemaxpar;
   fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);      }
   fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);      
   printf("begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);      fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
   fprintf(ficlog,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);      fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
          fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
   while((c=getc(ficpar))=='#' && c!= EOF){      
       while((c=getc(ficpar))=='#' && c!= EOF){
         ungetc(c,ficpar);
         fgets(line, MAXLINE, ficpar);
         puts(line);
         fputs(line,ficparo);
       }
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      
     puts(line);      fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
     fputs(line,ficparo);      fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
   }      fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
   ungetc(c,ficpar);      printf("begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
        fprintf(ficlog,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
       
   dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.;      while((c=getc(ficpar))=='#' && c!= EOF){
   dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.;        ungetc(c,ficpar);
         fgets(line, MAXLINE, ficpar);
   fscanf(ficpar,"pop_based=%d\n",&popbased);        puts(line);
   fprintf(ficparo,"pop_based=%d\n",popbased);           fputs(line,ficparo);
   fprintf(ficres,"pop_based=%d\n",popbased);         }
     
   while((c=getc(ficpar))=='#' && c!= EOF){  
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      
     puts(line);      
     fputs(line,ficparo);      dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.;
   }      dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.;
   ungetc(c,ficpar);      
       fscanf(ficpar,"pop_based=%d\n",&popbased);
   fscanf(ficpar,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj);      fprintf(ficparo,"pop_based=%d\n",popbased);   
   fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);      fprintf(ficres,"pop_based=%d\n",popbased);   
   printf("prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);      
   fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);      while((c=getc(ficpar))=='#' && c!= EOF){
   fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);        ungetc(c,ficpar);
   /* day and month of proj2 are not used but only year anproj2.*/        fgets(line, MAXLINE, ficpar);
         puts(line);
   while((c=getc(ficpar))=='#' && c!= EOF){        fputs(line,ficparo);
       }
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      
     puts(line);      fscanf(ficpar,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj);
     fputs(line,ficparo);      fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
   }      printf("prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
   ungetc(c,ficpar);      fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
       fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
   fscanf(ficpar,"popforecast=%d popfile=%s popfiledate=%lf/%lf/%lf last-popfiledate=%lf/%lf/%lf\n",&popforecast,popfile,&jpyram,&mpyram,&anpyram,&jpyram1,&mpyram1,&anpyram1);      /* day and month of proj2 are not used but only year anproj2.*/
   fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);      
   fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);      
       
   /*  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);*/      /*  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);*/
   /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/      /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
       
   replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */      replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */
   printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);      printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
       
   printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\      printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\
                model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\                   model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\
                jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);                   jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
          
   /*------------ free_vector  -------------*/     /*------------ free_vector  -------------*/
   /*  chdir(path); */     /*  chdir(path); */
     
   free_ivector(wav,1,imx);      free_ivector(wav,1,imx);
   free_imatrix(dh,1,lastpass-firstpass+1,1,imx);      free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
   free_imatrix(bh,1,lastpass-firstpass+1,1,imx);      free_imatrix(bh,1,lastpass-firstpass+1,1,imx);
   free_imatrix(mw,1,lastpass-firstpass+1,1,imx);         free_imatrix(mw,1,lastpass-firstpass+1,1,imx);   
   free_lvector(num,1,n);      free_lvector(num,1,n);
   free_vector(agedc,1,n);      free_vector(agedc,1,n);
   /*free_matrix(covar,0,NCOVMAX,1,n);*/      /*free_matrix(covar,0,NCOVMAX,1,n);*/
   /*free_matrix(covar,1,NCOVMAX,1,n);*/      /*free_matrix(covar,1,NCOVMAX,1,n);*/
   fclose(ficparo);      fclose(ficparo);
   fclose(ficres);      fclose(ficres);
   
   
   /*--------------- Prevalence limit  (stable prevalence) --------------*/      /*--------------- Prevalence limit  (stable prevalence) --------------*/
       
   strcpy(filerespl,"pl");      strcpy(filerespl,"pl");
   strcat(filerespl,fileres);      strcat(filerespl,fileres);
   if((ficrespl=fopen(filerespl,"w"))==NULL) {      if((ficrespl=fopen(filerespl,"w"))==NULL) {
     printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end;        printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end;
     fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end;        fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end;
   }      }
   printf("Computing stable prevalence: result on file '%s' \n", filerespl);      printf("Computing stable prevalence: result on file '%s' \n", filerespl);
   fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl);      fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl);
   fprintf(ficrespl,"#Stable prevalence \n");      fprintf(ficrespl,"#Stable prevalence \n");
   fprintf(ficrespl,"#Age ");      fprintf(ficrespl,"#Age ");
   for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);      for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
   fprintf(ficrespl,"\n");      fprintf(ficrespl,"\n");
       
   prlim=matrix(1,nlstate,1,nlstate);      prlim=matrix(1,nlstate,1,nlstate);
   
   agebase=ageminpar;      agebase=ageminpar;
   agelim=agemaxpar;      agelim=agemaxpar;
   ftolpl=1.e-10;      ftolpl=1.e-10;
   i1=cptcoveff;      i1=cptcoveff;
   if (cptcovn < 1){i1=1;}      if (cptcovn < 1){i1=1;}
   
   for(cptcov=1,k=0;cptcov<=i1;cptcov++){      for(cptcov=1,k=0;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){        for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
       k=k+1;          k=k+1;
       /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/          /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/
       fprintf(ficrespl,"\n#******");          fprintf(ficrespl,"\n#******");
       printf("\n#******");          printf("\n#******");
       fprintf(ficlog,"\n#******");          fprintf(ficlog,"\n#******");
       for(j=1;j<=cptcoveff;j++) {          for(j=1;j<=cptcoveff;j++) {
         fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);            fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);            printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);            fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
       }          }
       fprintf(ficrespl,"******\n");          fprintf(ficrespl,"******\n");
       printf("******\n");          printf("******\n");
       fprintf(ficlog,"******\n");          fprintf(ficlog,"******\n");
                   
       for (age=agebase; age<=agelim; age++){          for (age=agebase; age<=agelim; age++){
         prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);            prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
         fprintf(ficrespl,"%.0f ",age );            fprintf(ficrespl,"%.0f ",age );
         for(j=1;j<=cptcoveff;j++)            for(j=1;j<=cptcoveff;j++)
           fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);              fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         for(i=1; i<=nlstate;i++)            for(i=1; i<=nlstate;i++)
           fprintf(ficrespl," %.5f", prlim[i][i]);              fprintf(ficrespl," %.5f", prlim[i][i]);
         fprintf(ficrespl,"\n");            fprintf(ficrespl,"\n");
           }
       }        }
     }      }
   }      fclose(ficrespl);
   fclose(ficrespl);  
   
   /*------------- h Pij x at various ages ------------*/      /*------------- h Pij x at various ages ------------*/
       
   strcpy(filerespij,"pij");  strcat(filerespij,fileres);      strcpy(filerespij,"pij");  strcat(filerespij,fileres);
   if((ficrespij=fopen(filerespij,"w"))==NULL) {      if((ficrespij=fopen(filerespij,"w"))==NULL) {
     printf("Problem with Pij resultfile: %s\n", filerespij);goto end;        printf("Problem with Pij resultfile: %s\n", filerespij);goto end;
     fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij);goto end;        fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij);goto end;
   }      }
   printf("Computing pij: result on file '%s' \n", filerespij);      printf("Computing pij: result on file '%s' \n", filerespij);
   fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij);      fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij);
       
   stepsize=(int) (stepm+YEARM-1)/YEARM;      stepsize=(int) (stepm+YEARM-1)/YEARM;
   /*if (stepm<=24) stepsize=2;*/      /*if (stepm<=24) stepsize=2;*/
   
   agelim=AGESUP;      agelim=AGESUP;
   hstepm=stepsize*YEARM; /* Every year of age */      hstepm=stepsize*YEARM; /* Every year of age */
   hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */       hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ 
   
   /* hstepm=1;   aff par mois*/      /* hstepm=1;   aff par mois*/
   
   fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");      fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
   for(cptcov=1,k=0;cptcov<=i1;cptcov++){      for(cptcov=1,k=0;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){        for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
       k=k+1;          k=k+1;
       fprintf(ficrespij,"\n#****** ");          fprintf(ficrespij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++)           for(j=1;j<=cptcoveff;j++) 
         fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);            fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
       fprintf(ficrespij,"******\n");          fprintf(ficrespij,"******\n");
                   
       for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */          for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
         nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */             nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ 
         nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */            nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
   
         /*        nhstepm=nhstepm*YEARM; aff par mois*/            /*      nhstepm=nhstepm*YEARM; aff par mois*/
   
         p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);            p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
         oldm=oldms;savm=savms;            oldm=oldms;savm=savms;
         hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);              hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
         fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j=");            fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j=");
         for(i=1; i<=nlstate;i++)  
           for(j=1; j<=nlstate+ndeath;j++)  
             fprintf(ficrespij," %1d-%1d",i,j);  
         fprintf(ficrespij,"\n");  
         for (h=0; h<=nhstepm; h++){  
           fprintf(ficrespij,"%d %3.f %3.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );  
           for(i=1; i<=nlstate;i++)            for(i=1; i<=nlstate;i++)
             for(j=1; j<=nlstate+ndeath;j++)              for(j=1; j<=nlstate+ndeath;j++)
               fprintf(ficrespij," %.5f", p3mat[i][j][h]);                fprintf(ficrespij," %1d-%1d",i,j);
             fprintf(ficrespij,"\n");
             for (h=0; h<=nhstepm; h++){
               fprintf(ficrespij,"%d %3.f %3.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
               for(i=1; i<=nlstate;i++)
                 for(j=1; j<=nlstate+ndeath;j++)
                   fprintf(ficrespij," %.5f", p3mat[i][j][h]);
               fprintf(ficrespij,"\n");
             }
             free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           fprintf(ficrespij,"\n");            fprintf(ficrespij,"\n");
         }          }
         free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);  
         fprintf(ficrespij,"\n");  
       }        }
     }      }
   }  
   
   varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax);  
   
   fclose(ficrespij);      varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax);
   
   probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);      fclose(ficrespij);
   
   /*---------- Forecasting ------------------*/      probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
   /*if((stepm == 1) && (strcmp(model,".")==0)){*/      for(i=1;i<=AGESUP;i++)
   if(prevfcast==1){        for(j=1;j<=NCOVMAX;j++)
     /*    if(stepm ==1){*/          for(k=1;k<=NCOVMAX;k++)
             probs[i][j][k]=0.;
   
       /*---------- Forecasting ------------------*/
       /*if((stepm == 1) && (strcmp(model,".")==0)){*/
       if(prevfcast==1){
         /*    if(stepm ==1){*/
       prevforecast(fileres, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);        prevforecast(fileres, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
       /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/        /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/
 /*      }  */        /*      }  */
 /*      else{ */        /*      else{ */
 /*        erreur=108; */        /*        erreur=108; */
 /*        printf("Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */        /*        printf("Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
 /*        fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */        /*        fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
 /*      } */        /*      } */
   }      }
       
   
   /*---------- Health expectancies and variances ------------*/      /*---------- Health expectancies and variances ------------*/
   
   strcpy(filerest,"t");      strcpy(filerest,"t");
   strcat(filerest,fileres);      strcat(filerest,fileres);
   if((ficrest=fopen(filerest,"w"))==NULL) {      if((ficrest=fopen(filerest,"w"))==NULL) {
     printf("Problem with total LE resultfile: %s\n", filerest);goto end;        printf("Problem with total LE resultfile: %s\n", filerest);goto end;
     fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;        fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
   }      }
   printf("Computing Total LEs with variances: file '%s' \n", filerest);       printf("Computing Total LEs with variances: file '%s' \n", filerest); 
   fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest);       fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest); 
   
   
   strcpy(filerese,"e");      strcpy(filerese,"e");
   strcat(filerese,fileres);      strcat(filerese,fileres);
   if((ficreseij=fopen(filerese,"w"))==NULL) {      if((ficreseij=fopen(filerese,"w"))==NULL) {
     printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);        printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
     fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0);        fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
   }      }
   printf("Computing Health Expectancies: result on file '%s' \n", filerese);      printf("Computing Health Expectancies: result on file '%s' \n", filerese);
   fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);      fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
   
   strcpy(fileresv,"v");      strcpy(fileresv,"v");
   strcat(fileresv,fileres);      strcat(fileresv,fileres);
   if((ficresvij=fopen(fileresv,"w"))==NULL) {      if((ficresvij=fopen(fileresv,"w"))==NULL) {
     printf("Problem with variance resultfile: %s\n", fileresv);exit(0);        printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
     fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0);        fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0);
   }      }
   printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);      printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
   fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);      fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
   
   /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */      /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */
   prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);      prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
   /*  printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d,  mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\      /*  printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d,  mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\
 ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);          ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);
   */      */
   
   if (mobilav!=0) {      if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);        mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){        if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
       fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);          fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
       printf(" Error in movingaverage mobilav=%d\n",mobilav);          printf(" Error in movingaverage mobilav=%d\n",mobilav);
         }
     }      }
   }  
   
   for(cptcov=1,k=0;cptcov<=i1;cptcov++){      for(cptcov=1,k=0;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){        for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
       k=k+1;           k=k+1; 
       fprintf(ficrest,"\n#****** ");          fprintf(ficrest,"\n#****** ");
       for(j=1;j<=cptcoveff;j++)           for(j=1;j<=cptcoveff;j++) 
         fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);            fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
       fprintf(ficrest,"******\n");          fprintf(ficrest,"******\n");
   
       fprintf(ficreseij,"\n#****** ");          fprintf(ficreseij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++)           for(j=1;j<=cptcoveff;j++) 
         fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);            fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
       fprintf(ficreseij,"******\n");          fprintf(ficreseij,"******\n");
   
       fprintf(ficresvij,"\n#****** ");          fprintf(ficresvij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++)           for(j=1;j<=cptcoveff;j++) 
         fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);            fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
       fprintf(ficresvij,"******\n");          fprintf(ficresvij,"******\n");
   
       eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);          eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
       oldm=oldms;savm=savms;          oldm=oldms;savm=savms;
       evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov);            evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov);  
     
       vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);          vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
       oldm=oldms;savm=savms;          oldm=oldms;savm=savms;
       varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav);          varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav);
       if(popbased==1){          if(popbased==1){
         varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav);            varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav);
       }          }
   
     
       fprintf(ficrest,"#Total LEs with variances: e.. (std) ");          fprintf(ficrest,"#Total LEs with variances: e.. (std) ");
       for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);          for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
       fprintf(ficrest,"\n");          fprintf(ficrest,"\n");
   
       epj=vector(1,nlstate+1);          epj=vector(1,nlstate+1);
       for(age=bage; age <=fage ;age++){          for(age=bage; age <=fage ;age++){
         prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);            prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
         if (popbased==1) {            if (popbased==1) {
           if(mobilav ==0){              if(mobilav ==0){
             for(i=1; i<=nlstate;i++)                for(i=1; i<=nlstate;i++)
               prlim[i][i]=probs[(int)age][i][k];                  prlim[i][i]=probs[(int)age][i][k];
           }else{ /* mobilav */               }else{ /* mobilav */ 
             for(i=1; i<=nlstate;i++)                for(i=1; i<=nlstate;i++)
               prlim[i][i]=mobaverage[(int)age][i][k];                  prlim[i][i]=mobaverage[(int)age][i][k];
               }
           }            }
         }  
                   
         fprintf(ficrest," %4.0f",age);            fprintf(ficrest," %4.0f",age);
         for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){            for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
           for(i=1, epj[j]=0.;i <=nlstate;i++) {              for(i=1, epj[j]=0.;i <=nlstate;i++) {
             epj[j] += prlim[i][i]*eij[i][j][(int)age];                epj[j] += prlim[i][i]*eij[i][j][(int)age];
             /*  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/                /*  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
               }
               epj[nlstate+1] +=epj[j];
           }            }
           epj[nlstate+1] +=epj[j];  
         }  
   
         for(i=1, vepp=0.;i <=nlstate;i++)            for(i=1, vepp=0.;i <=nlstate;i++)
           for(j=1;j <=nlstate;j++)              for(j=1;j <=nlstate;j++)
             vepp += vareij[i][j][(int)age];                vepp += vareij[i][j][(int)age];
         fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));            fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
         for(j=1;j <=nlstate;j++){            for(j=1;j <=nlstate;j++){
           fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));              fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
             }
             fprintf(ficrest,"\n");
         }          }
         fprintf(ficrest,"\n");          free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
       }          free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
       free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);          free_vector(epj,1,nlstate+1);
       free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);        }
       free_vector(epj,1,nlstate+1);      }
     }      free_vector(weight,1,n);
   }      free_imatrix(Tvard,1,15,1,2);
   free_vector(weight,1,n);      free_imatrix(s,1,maxwav+1,1,n);
   free_imatrix(Tvard,1,15,1,2);      free_matrix(anint,1,maxwav,1,n); 
   free_imatrix(s,1,maxwav+1,1,n);      free_matrix(mint,1,maxwav,1,n);
   free_matrix(anint,1,maxwav,1,n);       free_ivector(cod,1,n);
   free_matrix(mint,1,maxwav,1,n);      free_ivector(tab,1,NCOVMAX);
   free_ivector(cod,1,n);      fclose(ficreseij);
   free_ivector(tab,1,NCOVMAX);      fclose(ficresvij);
   fclose(ficreseij);      fclose(ficrest);
   fclose(ficresvij);      fclose(ficpar);
   fclose(ficrest);    
   fclose(ficpar);      /*------- Variance of stable prevalence------*/   
     
   /*------- Variance of stable prevalence------*/         strcpy(fileresvpl,"vpl");
       strcat(fileresvpl,fileres);
   strcpy(fileresvpl,"vpl");      if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
   strcat(fileresvpl,fileres);        printf("Problem with variance of stable prevalence  resultfile: %s\n", fileresvpl);
   if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {        exit(0);
     printf("Problem with variance of stable prevalence  resultfile: %s\n", fileresvpl);      }
     exit(0);      printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl);
   }  
   printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl);      for(cptcov=1,k=0;cptcov<=i1;cptcov++){
         for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
   for(cptcov=1,k=0;cptcov<=i1;cptcov++){          k=k+1;
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){          fprintf(ficresvpl,"\n#****** ");
       k=k+1;          for(j=1;j<=cptcoveff;j++) 
       fprintf(ficresvpl,"\n#****** ");            fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
       for(j=1;j<=cptcoveff;j++)           fprintf(ficresvpl,"******\n");
         fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);  
       fprintf(ficresvpl,"******\n");  
               
       varpl=matrix(1,nlstate,(int) bage, (int) fage);          varpl=matrix(1,nlstate,(int) bage, (int) fage);
       oldm=oldms;savm=savms;          oldm=oldms;savm=savms;
       varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);          varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);
       free_matrix(varpl,1,nlstate,(int) bage, (int)fage);          free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
         }
     }      }
   }  
   
   fclose(ficresvpl);      fclose(ficresvpl);
   
   /*---------- End : free ----------------*/      /*---------- End : free ----------------*/
   free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);      if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);      free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
   free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);  
   free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);    }  /* mle==-3 arrives here for freeing */
         free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(covar,0,NCOVMAX,1,n);      free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(matcov,1,npar,1,npar);      free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
   /*free_vector(delti,1,npar);*/      free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);     
   free_matrix(agev,1,maxwav,1,imx);      free_matrix(covar,0,NCOVMAX,1,n);
   free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);      free_matrix(matcov,1,npar,1,npar);
   if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);      /*free_vector(delti,1,npar);*/
   free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);      free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); 
       free_matrix(agev,1,maxwav,1,imx);
       free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
   
       free_ivector(ncodemax,1,8);
       free_ivector(Tvar,1,15);
       free_ivector(Tprod,1,15);
       free_ivector(Tvaraff,1,15);
       free_ivector(Tage,1,15);
       free_ivector(Tcode,1,100);
   
   free_ivector(ncodemax,1,8);  
   free_ivector(Tvar,1,15);  
   free_ivector(Tprod,1,15);  
   free_ivector(Tvaraff,1,15);  
   free_ivector(Tage,1,15);  
   free_ivector(Tcode,1,100);  
   
   fflush(fichtm);    fflush(fichtm);
   fflush(ficgp);    fflush(ficgp);
Line 5192  ageminpar, agemax, s[lastpass][imx], age Line 5433  ageminpar, agemax, s[lastpass][imx], age
   /*------ End -----------*/    /*------ End -----------*/
   
   chdir(path);    chdir(path);
   strcpy(plotcmd,GNUPLOTPROGRAM);    strcpy(plotcmd,"\"");
     strcat(plotcmd,pathimach);
     strcat(plotcmd,GNUPLOTPROGRAM);
     strcat(plotcmd,"\"");
   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);
Line 5205  ageminpar, agemax, s[lastpass][imx], age Line 5449  ageminpar, agemax, s[lastpass][imx], age
     printf("\nType e to edit output files, g to graph 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"); */
     if (z[0] == 'e') system(optionfilehtm);      if (z[0] == 'e') {
         printf("Starting browser with: %s",optionfilehtm);fflush(stdout);
         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);
   }    }
Removed from v.1.95  
changed lines
  Added in v.1.102

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