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

Diff for /imach/src/imach.c between versions 1.10 and 1.27

version 1.10, 2001/05/09 14:25:42 version 1.27, 2002/02/28 17:49:07
Line 1 Line 1
      /* $Id$
 /*********************** Imach **************************************             Interpolated Markov Chain
   This program computes Healthy Life Expectancies from cross-longitudinal  
   data. Cross-longitudinal consist in a first survey ("cross") where    Short summary of the programme:
   individuals from different ages are interviewed on their health status   
   or degree of  disability. At least a second wave of interviews    This program computes Healthy Life Expectancies from
   ("longitudinal") should  measure each new individual health status.    cross-longitudinal data. Cross-longitudinal data consist in: -1- a
   Health expectancies are computed from the transistions observed between    first survey ("cross") where individuals from different ages are
   waves and are computed for each degree of severity of disability (number    interviewed on their health status or degree of disability (in the
   of life states). More degrees you consider, more time is necessary to    case of a health survey which is our main interest) -2- at least a
   reach the Maximum Likelihood of the parameters involved in the model.    second wave of interviews ("longitudinal") which measure each change
   The simplest model is the multinomial logistic model where pij is    (if any) in individual health status.  Health expectancies are
   the probabibility to be observed in state j at the second wave conditional    computed from the time spent in each health state according to a
   to be observed in state i at the first wave. Therefore the model is:    model. More health states you consider, more time is necessary to reach the
   log(pij/pii)= aij + bij*age+ cij*sex + etc , where 'age' is age and 'sex'    Maximum Likelihood of the parameters involved in the model.  The
   is a covariate. If you want to have a more complex model than "constant and    simplest model is the multinomial logistic model where pij is the
   age", you should modify the program where the markup    probabibility to be observed in state j at the second wave
     *Covariates have to be included here again* invites you to do it.    conditional to be observed in state i at the first wave. Therefore
   More covariates you add, less is the speed of the convergence.    the model is: log(pij/pii)= aij + bij*age+ cij*sex + etc , where
     'age' is age and 'sex' is a covariate. If you want to have a more
   The advantage that this computer programme claims, comes from that if the    complex model than "constant and age", you should modify the program
   delay between waves is not identical for each individual, or if some    where the markup *Covariates have to be included here again* invites
   individual missed an interview, the information is not rounded or lost, but    you to do it.  More covariates you add, slower the
   taken into account using an interpolation or extrapolation.    convergence.
   hPijx is the probability to be  
   observed in state i at age x+h conditional to the observed state i at age    The advantage of this computer programme, compared to a simple
   x. The delay 'h' can be split into an exact number (nh*stepm) of    multinomial logistic model, is clear when the delay between waves is not
   unobserved intermediate  states. This elementary transition (by month or    identical for each individual. Also, if a individual missed an
   quarter trimester, semester or year) is model as a multinomial logistic.    intermediate interview, the information is lost, but taken into
   The hPx matrix is simply the matrix product of nh*stepm elementary matrices    account using an interpolation or extrapolation.  
   and the contribution of each individual to the likelihood is simply hPijx.  
     hPijx is the probability to be observed in state i at age x+h
     conditional to the observed state i at age x. The delay 'h' can be
     split into an exact number (nh*stepm) of unobserved intermediate
     states. This elementary transition (by month or quarter trimester,
     semester or year) is model as a multinomial logistic.  The hPx
     matrix is simply the matrix product of nh*stepm elementary matrices
     and the contribution of each individual to the likelihood is simply
     hPijx.
   
   Also this programme outputs the covariance matrix of the parameters but also    Also this programme outputs the covariance matrix of the parameters but also
   of the life expectancies. It also computes the prevalence limits.    of the life expectancies. It also computes the prevalence limits.
Line 48 Line 56
 #include <unistd.h>  #include <unistd.h>
   
 #define MAXLINE 256  #define MAXLINE 256
   #define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"
 #define FILENAMELENGTH 80  #define FILENAMELENGTH 80
 /*#define DEBUG*/  /*#define DEBUG*/
 #define windows  #define windows
Line 67 Line 76
 #define AGEBASE 40  #define AGEBASE 40
   
   
   int erreur; /* Error number */
 int nvar;  int nvar;
 int cptcovn, cptcovage=0, cptcoveff=0,cptcov;  int cptcovn, cptcovage=0, cptcoveff=0,cptcov;
 int npar=NPARMAX;  int npar=NPARMAX;
 int nlstate=2; /* Number of live states */  int nlstate=2; /* Number of live states */
 int ndeath=1; /* Number of dead states */  int ndeath=1; /* Number of dead states */
 int ncovmodel, ncov;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */  int ncovmodel, ncov;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
   int popbased=0;
   
 int *wav; /* Number of waves for this individuual 0 is possible */  int *wav; /* Number of waves for this individuual 0 is possible */
 int maxwav; /* Maxim number of waves */  int maxwav; /* Maxim number of waves */
Line 83  int **dh; /* dh[mi][i] is number of step Line 94  int **dh; /* dh[mi][i] is number of step
 double jmean; /* Mean space between 2 waves */  double jmean; /* Mean space between 2 waves */
 double **oldm, **newm, **savm; /* Working pointers to matrices */  double **oldm, **newm, **savm; /* Working pointers to matrices */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */  double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
 FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest;  FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
 FILE *ficgp, *fichtm;  FILE *ficgp,*ficresprob,*ficpop;
 FILE *ficreseij;  FILE *ficreseij;
   char filerese[FILENAMELENGTH];    char filerese[FILENAMELENGTH];
  FILE  *ficresvij;   FILE  *ficresvij;
Line 112  FILE *ficreseij; Line 123  FILE *ficreseij;
 static double maxarg1,maxarg2;  static double maxarg1,maxarg2;
 #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)>(maxarg2)? (maxarg1):(maxarg2))  #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)>(maxarg2)? (maxarg1):(maxarg2))
 #define FMIN(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)<(maxarg2)? (maxarg1):(maxarg2))  #define FMIN(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)<(maxarg2)? (maxarg1):(maxarg2))
     
 #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a))  #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a))
 #define rint(a) floor(a+0.5)  #define rint(a) floor(a+0.5)
   
Line 127  int stepm; Line 138  int stepm;
 int m,nb;  int m,nb;
 int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;  int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;
 double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;  double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
 double **pmmij;  double **pmmij, ***probs, ***mobaverage;
   double dateintmean=0;
   
 double *weight;  double *weight;
 int **s; /* Status */  int **s; /* Status */
Line 138  double ftol=FTOL; /* Tolerance for compu Line 150  double ftol=FTOL; /* Tolerance for compu
 double ftolhess; /* Tolerance for computing hessian */  double ftolhess; /* Tolerance for computing hessian */
   
 /**************** split *************************/  /**************** split *************************/
 static  int split( char *path, char *dirc, char *name )  static  int split( char *path, char *dirc, char *name, char *ext, char *finame )
 {  {
    char *s;                             /* pointer */     char *s;                             /* pointer */
    int  l1, l2;                         /* length counters */     int  l1, l2;                         /* length counters */
   
    l1 = strlen( path );                 /* length of path */     l1 = strlen( path );                 /* length of path */
    if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );     if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
   #ifdef windows
    s = strrchr( path, '\\' );           /* find last / */     s = strrchr( path, '\\' );           /* find last / */
   #else
      s = strrchr( path, '/' );            /* find last / */
   #endif
    if ( s == NULL ) {                   /* no directory, so use current */     if ( s == NULL ) {                   /* no directory, so use current */
 #if     defined(__bsd__)                /* get current working directory */  #if     defined(__bsd__)                /* get current working directory */
       extern char       *getwd( );        extern char       *getwd( );
Line 168  static int split( char *path, char *dirc Line 184  static int split( char *path, char *dirc
       dirc[l1-l2] = 0;                  /* add zero */        dirc[l1-l2] = 0;                  /* add zero */
    }     }
    l1 = strlen( dirc );                 /* length of directory */     l1 = strlen( dirc );                 /* length of directory */
   #ifdef windows
    if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }     if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
   #else
      if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }
   #endif
      s = strrchr( name, '.' );            /* find last / */
      s++;
      strcpy(ext,s);                       /* save extension */
      l1= strlen( name);
      l2= strlen( s)+1;
      strncpy( finame, name, l1-l2);
      finame[l1-l2]= 0;
    return( 0 );                         /* we're done */     return( 0 );                         /* we're done */
 }  }
   
Line 693  double **prevalim(double **prlim, int nl Line 720  double **prevalim(double **prlim, int nl
   }    }
 }  }
   
 /*************** transition probabilities **********/  /*************** transition probabilities ***************/
   
 double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )  double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
 {  {
Line 719  double **pmij(double **ps, double *cov, Line 746  double **pmij(double **ps, double *cov,
       ps[i][j]=s2;        ps[i][j]=s2;
     }      }
   }    }
       /*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++)
Line 740  double **pmij(double **ps, double *cov, Line 769  double **pmij(double **ps, double *cov,
     }      }
   }    }
   
   
   /*   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("%lf ",ps[ii][jj]);
Line 757  double **pmij(double **ps, double *cov, Line 787  double **pmij(double **ps, double *cov,
   
 double **matprod2(double **out, double **in,long nrl, long nrh, long ncl, long nch, long ncolol, long ncoloh, double **b)  double **matprod2(double **out, double **in,long nrl, long nrh, long ncl, long nch, long ncolol, long ncoloh, double **b)
 {  {
   /* Computes the matric product of in(1,nrh-nrl+1)(1,nch-ncl+1) times    /* Computes the matrix product of in(1,nrh-nrl+1)(1,nch-ncl+1) times
      b(1,nch-ncl+1)(1,ncoloh-ncolol+1) into out(...) */       b(1,nch-ncl+1)(1,ncoloh-ncolol+1) into out(...) */
   /* in, b, out are matrice of pointers which should have been initialized    /* in, b, out are matrice of pointers which should have been initialized
      before: only the contents of out is modified. The function returns       before: only the contents of out is modified. The function returns
Line 804  double ***hpxij(double ***po, int nhstep Line 834  double ***hpxij(double ***po, int nhstep
       cov[1]=1.;        cov[1]=1.;
       cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;        cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
       for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];        for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
 for (k=1; k<=cptcovage;k++)        for (k=1; k<=cptcovage;k++)
         cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];          cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
    for (k=1; k<=cptcovprod;k++)        for (k=1; k<=cptcovprod;k++)
         cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];          cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
   
   
Line 896  void mlikeli(FILE *ficres,double p[], in Line 926  void mlikeli(FILE *ficres,double p[], in
   powell(p,xi,npar,ftol,&iter,&fret,func);    powell(p,xi,npar,ftol,&iter,&fret,func);
   
    printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));     printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
   fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f ",iter,func(p));    fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
   
 }  }
   
Line 913  void hesscov(double **matcov, double p[] Line 943  void hesscov(double **matcov, double p[]
   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) ;
   
   
   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 921  void hesscov(double **matcov, double p[] Line 950  void hesscov(double **matcov, double p[]
     printf("%d",i);fflush(stdout);      printf("%d",i);fflush(stdout);
     hess[i][i]=hessii(p,ftolhess,i,delti);      hess[i][i]=hessii(p,ftolhess,i,delti);
     /*printf(" %f ",p[i]);*/      /*printf(" %f ",p[i]);*/
       /*printf(" %lf ",hess[i][i]);*/
   }    }
    
   for (i=1;i<=npar;i++) {    for (i=1;i<=npar;i++) {
     for (j=1;j<=npar;j++)  {      for (j=1;j<=npar;j++)  {
       if (j>i) {        if (j>i) {
         printf(".%d%d",i,j);fflush(stdout);          printf(".%d%d",i,j);fflush(stdout);
         hess[i][j]=hessij(p,delti,i,j);          hess[i][j]=hessij(p,delti,i,j);
         hess[j][i]=hess[i][j];          hess[j][i]=hess[i][j];    
           /*printf(" %lf ",hess[i][j]);*/
       }        }
     }      }
   }    }
Line 1032  double hessii( double x[], double delta, Line 1063  double hessii( double x[], double delta,
     }      }
   }    }
   delti[theta]=delts;    delti[theta]=delts;
   return res;    return res;
     
 }  }
   
Line 1145  void lubksb(double **a, int n, int *indx Line 1176  void lubksb(double **a, int n, int *indx
 }  }
   
 /************ Frequencies ********************/  /************ Frequencies ********************/
 void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax)  void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2)
 {  /* Some frequencies */  {  /* Some frequencies */
     
   int i, m, jk, k1, i1, j1, bool, z1,z2,j;    int i, m, jk, k1,i1, j1, bool, z1,z2,j;
   double ***freq; /* Frequencies */    double ***freq; /* Frequencies */
   double *pp;    double *pp;
   double pos;    double pos, k2, dateintsum=0,k2cpt=0;
   FILE *ficresp;    FILE *ficresp;
   char fileresp[FILENAMELENGTH];    char fileresp[FILENAMELENGTH];
   
   pp=vector(1,nlstate);    pp=vector(1,nlstate);
     probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
   strcpy(fileresp,"p");    strcpy(fileresp,"p");
   strcat(fileresp,fileres);    strcat(fileresp,fileres);
   if((ficresp=fopen(fileresp,"w"))==NULL) {    if((ficresp=fopen(fileresp,"w"))==NULL) {
Line 1178  void  freqsummary(char fileres[], int ag Line 1209  void  freqsummary(char fileres[], int ag
          for (jk=-1; jk<=nlstate+ndeath; jk++)             for (jk=-1; jk<=nlstate+ndeath; jk++)  
            for(m=agemin; m <= agemax+3; m++)             for(m=agemin; m <= agemax+3; m++)
              freq[i][jk][m]=0;               freq[i][jk][m]=0;
          
           dateintsum=0;
           k2cpt=0;
        for (i=1; i<=imx; i++) {         for (i=1; i<=imx; i++) {
          bool=1;           bool=1;
          if  (cptcovn>0) {           if  (cptcovn>0) {
Line 1186  void  freqsummary(char fileres[], int ag Line 1219  void  freqsummary(char fileres[], int ag
              if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])               if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
                bool=0;                 bool=0;
          }           }
           if (bool==1) {           if (bool==1) {
            for(m=firstpass; m<=lastpass-1; m++){             for(m=firstpass; m<=lastpass; m++){
              if(agev[m][i]==0) agev[m][i]=agemax+1;               k2=anint[m][i]+(mint[m][i]/12.);
              if(agev[m][i]==1) agev[m][i]=agemax+2;               if ((k2>=dateprev1) && (k2<=dateprev2)) {
              freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];                 if(agev[m][i]==0) agev[m][i]=agemax+1;
              freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];                 if(agev[m][i]==1) agev[m][i]=agemax+2;
                  freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
                  freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];
                  if ((agev[m][i]>1) && (agev[m][i]< (agemax+3))) {
                    dateintsum=dateintsum+k2;
                    k2cpt++;
                  }
   
                }
            }             }
          }           }
        }         }
          
          fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
   
         if  (cptcovn>0) {          if  (cptcovn>0) {
          fprintf(ficresp, "\n#********** Variable ");           fprintf(ficresp, "\n#********** Variable ");
          for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);           for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
Line 1211  void  freqsummary(char fileres[], int ag Line 1255  void  freqsummary(char fileres[], int ag
       printf("Age %d", i);        printf("Age %d", i);
     for(jk=1; jk <=nlstate ; jk++){      for(jk=1; jk <=nlstate ; jk++){
       for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)        for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
         pp[jk] += freq[jk][m][i];          pp[jk] += freq[jk][m][i];
     }      }
     for(jk=1; jk <=nlstate ; jk++){      for(jk=1; jk <=nlstate ; jk++){
       for(m=-1, pos=0; m <=0 ; m++)        for(m=-1, pos=0; m <=0 ; m++)
Line 1221  void  freqsummary(char fileres[], int ag Line 1265  void  freqsummary(char fileres[], int ag
       else        else
         printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);          printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
     }      }
     for(jk=1; jk <=nlstate ; jk++){  
       for(m=1, pp[jk]=0; m <=nlstate+ndeath; m++)       for(jk=1; jk <=nlstate ; jk++){
         for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
         pp[jk] += freq[jk][m][i];          pp[jk] += freq[jk][m][i];
     }       }
   
     for(jk=1,pos=0; jk <=nlstate ; jk++)      for(jk=1,pos=0; jk <=nlstate ; jk++)
       pos += pp[jk];        pos += pp[jk];
     for(jk=1; jk <=nlstate ; jk++){      for(jk=1; jk <=nlstate ; jk++){
Line 1233  void  freqsummary(char fileres[], int ag Line 1279  void  freqsummary(char fileres[], int ag
       else        else
         printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);          printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
       if( i <= (int) agemax){        if( i <= (int) agemax){
         if(pos>=1.e-5)          if(pos>=1.e-5){
           fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);            fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);
             probs[i][jk][j1]= pp[jk]/pos;
             /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/
           }
       else        else
           fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);            fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);
       }        }
Line 1248  void  freqsummary(char fileres[], int ag Line 1297  void  freqsummary(char fileres[], int ag
     }      }
     }      }
  }   }
     dateintmean=dateintsum/k2cpt;
     
   fclose(ficresp);    fclose(ficresp);
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);    free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);    free_vector(pp,1,nlstate);
   
     /* End of Freq */
   }
   
   /************ Prevalence ********************/
   void prevalence(int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, double calagedate)
   {  /* Some frequencies */
    
     int i, m, jk, k1, i1, j1, bool, z1,z2,j;
     double ***freq; /* Frequencies */
     double *pp;
     double pos, k2;
   
     pp=vector(1,nlstate);
     probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
    
     freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
     j1=0;
    
     j=cptcoveff;
     if (cptcovn<1) {j=1;ncodemax[1]=1;}
    
    for(k1=1; k1<=j;k1++){
       for(i1=1; i1<=ncodemax[k1];i1++){
         j1++;
    
         for (i=-1; i<=nlstate+ndeath; i++)  
           for (jk=-1; jk<=nlstate+ndeath; jk++)  
             for(m=agemin; m <= agemax+3; m++)
               freq[i][jk][m]=0;
        
         for (i=1; i<=imx; i++) {
           bool=1;
           if  (cptcovn>0) {
             for (z1=1; z1<=cptcoveff; z1++)
               if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
                 bool=0;
           }
           if (bool==1) {
             for(m=firstpass; m<=lastpass; m++){
               k2=anint[m][i]+(mint[m][i]/12.);
               if ((k2>=dateprev1) && (k2<=dateprev2)) {
                 if(agev[m][i]==0) agev[m][i]=agemax+1;
                 if(agev[m][i]==1) agev[m][i]=agemax+2;
                 freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
                 freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i];  
               }
             }
           }
         }
        
           for(i=(int)agemin; i <= (int)agemax+3; i++){
             for(jk=1; jk <=nlstate ; jk++){
               for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
                 pp[jk] += freq[jk][m][i];
             }
             for(jk=1; jk <=nlstate ; jk++){
               for(m=-1, pos=0; m <=0 ; m++)
               pos += freq[jk][m][i];
           }
          
            for(jk=1; jk <=nlstate ; jk++){
              for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
                pp[jk] += freq[jk][m][i];
            }
            
            for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];
   
            for(jk=1; jk <=nlstate ; jk++){          
              if( i <= (int) agemax){
                if(pos>=1.e-5){
                  probs[i][jk][j1]= pp[jk]/pos;
                }
              }
            }
            
           }
       }
     }
    
    
     free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
     free_vector(pp,1,nlstate);
    
 }  /* End of Freq */  }  /* End of Freq */
   
 /************* Waves Concatenation ***************/  /************* Waves Concatenation ***************/
Line 1270  void  concatwav(int wav[], int **dh, int Line 1403  void  concatwav(int wav[], int **dh, int
   /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;    /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
      double sum=0., jmean=0.;*/       double sum=0., jmean=0.;*/
   
 int j, k=0,jk, ju, jl;    int j, k=0,jk, ju, jl;
      double sum=0.;    double sum=0.;
 jmin=1e+5;    jmin=1e+5;
  jmax=-1;    jmax=-1;
 jmean=0.;    jmean=0.;
   for(i=1; i<=imx; i++){    for(i=1; i<=imx; i++){
     mi=0;      mi=0;
     m=firstpass;      m=firstpass;
Line 1309  jmean=0.; Line 1442  jmean=0.;
           if(j==0) j=1;  /* Survives at least one month after exam */            if(j==0) j=1;  /* Survives at least one month after exam */
           k=k+1;            k=k+1;
           if (j >= jmax) jmax=j;            if (j >= jmax) jmax=j;
           else if (j <= jmin)jmin=j;            if (j <= jmin) jmin=j;
           sum=sum+j;            sum=sum+j;
           if (j<0) printf("j=%d num=%d ",j,i);            /* if (j<10) printf("j=%d num=%d ",j,i); */
           }            }
         }          }
         else{          else{
Line 1319  jmean=0.; Line 1452  jmean=0.;
           k=k+1;            k=k+1;
           if (j >= jmax) jmax=j;            if (j >= jmax) jmax=j;
           else if (j <= jmin)jmin=j;            else if (j <= jmin)jmin=j;
             /*   if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
           sum=sum+j;            sum=sum+j;
         }          }
         jk= j/stepm;          jk= j/stepm;
Line 1335  jmean=0.; Line 1469  jmean=0.;
   }    }
   jmean=sum/k;    jmean=sum/k;
   printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);    printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
 }   }
 /*********** Tricode ****************************/  /*********** Tricode ****************************/
 void tricode(int *Tvar, int **nbcode, int imx)  void tricode(int *Tvar, int **nbcode, int imx)
 {  {
Line 1373  void tricode(int *Tvar, int **nbcode, in Line 1507  void tricode(int *Tvar, int **nbcode, in
   
  for (k=0; k<19; k++) Ndum[k]=0;   for (k=0; k<19; k++) Ndum[k]=0;
   
  for (i=1; i<=ncovmodel; i++) {   for (i=1; i<=ncovmodel-2; i++) {
       ij=Tvar[i];        ij=Tvar[i];
       Ndum[ij]++;        Ndum[ij]++;
     }      }
Line 1448  void varevsij(char fileres[], double *** Line 1582  void varevsij(char fileres[], double ***
   double **dnewm,**doldm;    double **dnewm,**doldm;
   int i, j, nhstepm, hstepm, h;    int i, j, nhstepm, hstepm, h;
   int k, cptcode;    int k, cptcode;
    double *xp;    double *xp;
   double **gp, **gm;    double **gp, **gm;
   double ***gradg, ***trgradg;    double ***gradg, ***trgradg;
   double ***p3mat;    double ***p3mat;
Line 1484  void varevsij(char fileres[], double *** Line 1618  void varevsij(char fileres[], double ***
       }        }
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);          hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
   
         if (popbased==1) {
           for(i=1; i<=nlstate;i++)
             prlim[i][i]=probs[(int)age][i][ij];
         }
    
       for(j=1; j<= nlstate; j++){        for(j=1; j<= nlstate; j++){
         for(h=0; h<=nhstepm; h++){          for(h=0; h<=nhstepm; h++){
           for(i=1, gp[h][j]=0.;i<=nlstate;i++)            for(i=1, gp[h][j]=0.;i<=nlstate;i++)
Line 1495  void varevsij(char fileres[], double *** Line 1635  void varevsij(char fileres[], double ***
         xp[i] = x[i] - (i==theta ?delti[theta]:0);          xp[i] = x[i] - (i==theta ?delti[theta]:0);
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);          hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
    
         if (popbased==1) {
           for(i=1; i<=nlstate;i++)
             prlim[i][i]=probs[(int)age][i][ij];
         }
   
       for(j=1; j<= nlstate; j++){        for(j=1; j<= nlstate; j++){
         for(h=0; h<=nhstepm; h++){          for(h=0; h<=nhstepm; h++){
           for(i=1, gm[h][j]=0.;i<=nlstate;i++)            for(i=1, gm[h][j]=0.;i<=nlstate;i++)
             gm[h][j] += prlim[i][i]*p3mat[i][j][h];              gm[h][j] += prlim[i][i]*p3mat[i][j][h];
         }          }
       }        }
   
       for(j=1; j<= nlstate; j++)        for(j=1; j<= nlstate; j++)
         for(h=0; h<=nhstepm; h++){          for(h=0; h<=nhstepm; h++){
           gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];            gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
Line 1540  void varevsij(char fileres[], double *** Line 1687  void varevsij(char fileres[], double ***
     free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);      free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   } /* End age */    } /* End age */
     
   free_vector(xp,1,npar);    free_vector(xp,1,npar);
   free_matrix(doldm,1,nlstate,1,npar);    free_matrix(doldm,1,nlstate,1,npar);
   free_matrix(dnewm,1,nlstate,1,nlstate);    free_matrix(dnewm,1,nlstate,1,nlstate);
Line 1630  void varprevlim(char fileres[], double * Line 1777  void varprevlim(char fileres[], double *
   
 }  }
   
   /************ Variance of one-step probabilities  ******************/
   void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij)
 /***********************************************/  
 /**************** Main Program *****************/  
 /***********************************************/  
   
 /*int main(int argc, char *argv[])*/  
 int main()  
 {  {
     int i, j;
     int k=0, cptcode;
     double **dnewm,**doldm;
     double *xp;
     double *gp, *gm;
     double **gradg, **trgradg;
     double age,agelim, cov[NCOVMAX];
     int theta;
     char fileresprob[FILENAMELENGTH];
   
   int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;    strcpy(fileresprob,"prob");
   double agedeb, agefin,hf;    strcat(fileresprob,fileres);
   double agemin=1.e20, agemax=-1.e20;    if((ficresprob=fopen(fileresprob,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", fileresprob);
   double fret;    }
   double **xi,tmp,delta;    printf("Computing variance of one-step probabilities: result on file '%s' \n",fileresprob);
   
   double dum; /* Dummy variable */  
   double ***p3mat;  
   int *indx;  
   char line[MAXLINE], linepar[MAXLINE];  
   char title[MAXLINE];  
   char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];  
   char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];  
   char filerest[FILENAMELENGTH];  
   char fileregp[FILENAMELENGTH];  
   char path[80],pathc[80],pathcd[80],pathtot[80],model[20];  
   int firstobs=1, lastobs=10;  
   int sdeb, sfin; /* Status at beginning and end */  
   int c,  h , cpt,l;  
   int ju,jl, mi;  
   int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;  
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;  
     
   int hstepm, nhstepm;  
   double bage, fage, age, agelim, agebase;  
   double ftolpl=FTOL;  
   double **prlim;  
   double *severity;  
   double ***param; /* Matrix of parameters */  
   double  *p;  
   double **matcov; /* Matrix of covariance */  
   double ***delti3; /* Scale */  
   double *delti; /* Scale */  
   double ***eij, ***vareij;  
   double **varpl; /* Variances of prevalence limits by age */  
   double *epj, vepp;  
   char version[80]="Imach version 64b, May 2001, INED-EUROREVES ";  
   char *alph[]={"a","a","b","c","d","e"}, str[4];  
   
   char z[1]="c", occ;    xp=vector(1,npar);
 #include <sys/time.h>    dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
 #include <time.h>    doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));
   char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];  
   /* long total_usecs;  
   struct timeval start_time, end_time;  
     
   gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */    cov[1]=1;
     for (age=bage; age<=fage; age ++){
       cov[2]=age;
       gradg=matrix(1,npar,1,9);
       trgradg=matrix(1,9,1,npar);
       gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
       gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
      
       for(theta=1; theta <=npar; theta++){
         for(i=1; i<=npar; i++)
           xp[i] = x[i] + (i==theta ?delti[theta]:0);
        
         pmij(pmmij,cov,ncovmodel,xp,nlstate);
      
         k=0;
         for(i=1; i<= (nlstate+ndeath); i++){
           for(j=1; j<=(nlstate+ndeath);j++){
              k=k+1;
             gp[k]=pmmij[i][j];
           }
         }
   
         for(i=1; i<=npar; i++)
           xp[i] = x[i] - (i==theta ?delti[theta]:0);
      
   
   printf("\nIMACH, Version 0.64b");        pmij(pmmij,cov,ncovmodel,xp,nlstate);
   printf("\nEnter the parameter file name: ");        k=0;
         for(i=1; i<=(nlstate+ndeath); i++){
           for(j=1; j<=(nlstate+ndeath);j++){
             k=k+1;
             gm[k]=pmmij[i][j];
           }
         }
        
          for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++)
              gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];  
       }
   
 #ifdef windows       for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)
   scanf("%s",pathtot);        for(theta=1; theta <=npar; theta++)
   getcwd(pathcd, size);        trgradg[j][theta]=gradg[theta][j];
   /*cygwin_split_path(pathtot,path,optionfile);   
     printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/       matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);
   /* cutv(path,optionfile,pathtot,'\\');*/       matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);
   
 split(pathtot, path,optionfile);       pmij(pmmij,cov,ncovmodel,x,nlstate);
   chdir(path);  
   replace(pathc,path);  
 #endif  
 #ifdef unix  
   scanf("%s",optionfile);  
 #endif  
   
 /*-------- arguments in the command line --------*/       k=0;
        for(i=1; i<=(nlstate+ndeath); i++){
          for(j=1; j<=(nlstate+ndeath);j++){
            k=k+1;
            gm[k]=pmmij[i][j];
           }
        }
        
        /*printf("\n%d ",(int)age);
        for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
          
   
   strcpy(fileres,"r");         printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
   strcat(fileres, optionfile);       }*/
   
   /*---------arguments file --------*/    fprintf(ficresprob,"\n%d ",(int)age);
   
   if((ficpar=fopen(optionfile,"r"))==NULL)    {    for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
     printf("Problem with optionfile %s\n",optionfile);      if (i== 2) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);
     goto end;  if (i== 4) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);
   }    }
   
   strcpy(filereso,"o");      free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
   strcat(filereso,fileres);      free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
   if((ficparo=fopen(filereso,"w"))==NULL) {      free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
     printf("Problem with Output resultfile: %s\n", filereso);goto end;      free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
   }  }
    free_vector(xp,1,npar);
   fclose(ficresprob);
   
   /* Reads comments: lines beginning with '#' */  }
   while((c=getc(ficpar))=='#' && c!= EOF){  
     ungetc(c,ficpar);  /******************* Printing html file ***********/
     fgets(line, MAXLINE, ficpar);  void printinghtml(char fileres[], char title[], char datafile[], int firstpass, int lastpass, int stepm, int weightopt, char model[],int imx,int jmin, int jmax, double jmeanint,char optionfile[],char optionfilehtm[] ){
     puts(line);    int jj1, k1, i1, cpt;
     fputs(line,ficparo);    FILE *fichtm;
     /*char optionfilehtm[FILENAMELENGTH];*/
   
     strcpy(optionfilehtm,optionfile);
     strcat(optionfilehtm,".htm");
     if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {
       printf("Problem with %s \n",optionfilehtm), exit(0);
   }    }
   ungetc(c,ficpar);  
   
   fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);   fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.7 </font> <hr size=\"2\" color=\"#EC5E5E\">
   printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt,model);  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>
   fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt,model);  
   
   covar=matrix(0,NCOVMAX,1,n);  Total number of observations=%d <br>
   cptcovn=0;  Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>
   if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;  <hr  size=\"2\" color=\"#EC5E5E\">
   <li>Outputs files<br><br>\n
           - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n
   - Estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>
           - Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>
           - Transition probabilities: <a href=\"pij%s\">pij%s</a><br>
           - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>
           - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>
           - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>
           - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>
           - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>
           - Prevalences and population forecasting: <a href=\"f%s\">f%s</a> <br>
           <br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
    
   fprintf(fichtm," <li>Graphs</li><p>");
   
   ncovmodel=2+cptcovn;   m=cptcoveff;
   nvar=ncovmodel-1; /* Suppressing age as a basic covariate */   if (cptcovn < 1) {m=1;ncodemax[1]=1;}
    
   /* Read guess parameters */  
   /* Reads comments: lines beginning with '#' */  
   while((c=getc(ficpar))=='#' && c!= EOF){  
     ungetc(c,ficpar);  
     fgets(line, MAXLINE, ficpar);  
     puts(line);  
     fputs(line,ficparo);  
   }  
   ungetc(c,ficpar);  
    
   param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);  
     for(i=1; i <=nlstate; i++)  
     for(j=1; j <=nlstate+ndeath-1; j++){  
       fscanf(ficpar,"%1d%1d",&i1,&j1);  
       fprintf(ficparo,"%1d%1d",i1,j1);  
       printf("%1d%1d",i,j);  
       for(k=1; k<=ncovmodel;k++){  
         fscanf(ficpar," %lf",&param[i][j][k]);  
         printf(" %lf",param[i][j][k]);  
         fprintf(ficparo," %lf",param[i][j][k]);  
       }  
       fscanf(ficpar,"\n");  
       printf("\n");  
       fprintf(ficparo,"\n");  
     }  
    
   npar= (nlstate+ndeath-1)*nlstate*ncovmodel;  
   p=param[1][1];  
    
   /* Reads comments: lines beginning with '#' */  
   while((c=getc(ficpar))=='#' && c!= EOF){  
     ungetc(c,ficpar);  
     fgets(line, MAXLINE, ficpar);  
     puts(line);  
     fputs(line,ficparo);  
   }  
   ungetc(c,ficpar);  
   
   delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);   jj1=0;
   delti=vector(1,npar); /* Scale of each paramater (output from hesscov) */   for(k1=1; k1<=m;k1++){
   for(i=1; i <=nlstate; i++){     for(i1=1; i1<=ncodemax[k1];i1++){
     for(j=1; j <=nlstate+ndeath-1; j++){         jj1++;
          if (cptcovn > 0) {
            fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
            for (cpt=1; cpt<=cptcoveff;cpt++)
              fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
            fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
          }
          fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>
   <img src=\"pe%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);    
          for(cpt=1; cpt<nlstate;cpt++){
            fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.gif<br>
   <img src=\"p%s%d%d.gif\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
          }
       for(cpt=1; cpt<=nlstate;cpt++) {
          fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
   interval) in state (%d): v%s%d%d.gif <br>
   <img src=\"v%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  
        }
        for(cpt=1; cpt<=nlstate;cpt++) {
           fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.gif <br>
   <img src=\"exp%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
        }
        fprintf(fichtm,"\n<br>- Total life expectancy by age and
   health expectancies in states (1) and (2): e%s%d.gif<br>
   <img src=\"e%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
   fprintf(fichtm,"\n</body>");
      }
      }
   fclose(fichtm);
   }
   
   /******************* Gnuplot file **************/
   void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double agemin, double agemax, double fage , char pathc[], double p[]){
   
     int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
   
     strcpy(optionfilegnuplot,optionfilefiname);
     strcat(optionfilegnuplot,".plt");
     if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
       printf("Problem with file %s",optionfilegnuplot);
     }
   
   #ifdef windows
       fprintf(ficgp,"cd \"%s\" \n",pathc);
   #endif
   m=pow(2,cptcoveff);
    
    /* 1eme*/
     for (cpt=1; cpt<= nlstate ; cpt ++) {
      for (k1=1; k1<= m ; k1 ++) {
   
   #ifdef windows
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",agemin,fage,fileres,k1-1,k1-1);
   #endif
   #ifdef unix
   fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",agemin,fage,fileres);
   #endif
   
   for (i=1; i<= nlstate ; i ++) {
     if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
     else fprintf(ficgp," \%%*lf (\%%*lf)");
   }
       fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1);
       for (i=1; i<= nlstate ; i ++) {
     if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
     else fprintf(ficgp," \%%*lf (\%%*lf)");
   }
     fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"vpl%s\" every :::%d::%d u 1:($2-2*$3) \"\%%lf",fileres,k1-1,k1-1);
        for (i=1; i<= nlstate ; i ++) {
     if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
     else fprintf(ficgp," \%%*lf (\%%*lf)");
   }  
        fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));
   #ifdef unix
   fprintf(ficgp,"\nset ter gif small size 400,300");
   #endif
   fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
      }
     }
     /*2 eme*/
   
     for (k1=1; k1<= m ; k1 ++) {
       fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",agemin,fage);
      
       for (i=1; i<= nlstate+1 ; i ++) {
         k=2*i;
         fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:2 \"\%%lf",fileres,k1-1,k1-1);
         for (j=1; j<= nlstate+1 ; j ++) {
     if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
     else fprintf(ficgp," \%%*lf (\%%*lf)");
   }  
         if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");
         else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);
       fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",fileres,k1-1,k1-1);
         for (j=1; j<= nlstate+1 ; j ++) {
           if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
           else fprintf(ficgp," \%%*lf (\%%*lf)");
   }  
         fprintf(ficgp,"\" t\"\" w l 0,");
        fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",fileres,k1-1,k1-1);
         for (j=1; j<= nlstate+1 ; j ++) {
     if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
     else fprintf(ficgp," \%%*lf (\%%*lf)");
   }  
         if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");
         else fprintf(ficgp,"\" t\"\" w l 0,");
       }
       fprintf(ficgp,"\nset out \"e%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),k1);
     }
    
     /*3eme*/
   
     for (k1=1; k1<= m ; k1 ++) {
       for (cpt=1; cpt<= nlstate ; cpt ++) {
         k=2+nlstate*(cpt-1);
         fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",agemin,fage,fileres,k1-1,k1-1,k,cpt);
         for (i=1; i< nlstate ; i ++) {
           fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+i,cpt,i+1);
         }
         fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
       }
       }
    
     /* CV preval stat */
       for (k1=1; k1<= m ; k1 ++) {
       for (cpt=1; cpt<nlstate ; cpt ++) {
         k=3;
         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",agemin,agemax,fileres,k1,k+cpt+1,k+1);
   
         for (i=1; i< nlstate ; i ++)
           fprintf(ficgp,"+$%d",k+i+1);
         fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);
        
         l=3+(nlstate+ndeath)*cpt;
         fprintf(ficgp,",\"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",fileres,k1,l+cpt+1,l+1);
         for (i=1; i< nlstate ; i ++) {
           l=3+(nlstate+ndeath)*cpt;
           fprintf(ficgp,"+$%d",l+i+1);
         }
         fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);  
         fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
       }
     }  
    
     /* proba elementaires */
      for(i=1,jk=1; i <=nlstate; i++){
       for(k=1; k <=(nlstate+ndeath); k++){
         if (k != i) {
           for(j=1; j <=ncovmodel; j++){
          
             fprintf(ficgp,"p%d=%f ",jk,p[jk]);
             jk++;
             fprintf(ficgp,"\n");
           }
         }
       }
       }
   
       for(jk=1; jk <=m; jk++) {
     fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",agemin,agemax);
      i=1;
      for(k2=1; k2<=nlstate; k2++) {
        k3=i;
        for(k=1; k<=(nlstate+ndeath); k++) {
          if (k != k2){
           fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
   ij=1;
           for(j=3; j <=ncovmodel; j++) {
             if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
               fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
               ij++;
             }
             else
             fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
           }
             fprintf(ficgp,")/(1");
          
           for(k1=1; k1 <=nlstate; k1++){  
             fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
   ij=1;
             for(j=3; j <=ncovmodel; j++){
             if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
               fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
               ij++;
             }
             else
               fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
             }
             fprintf(ficgp,")");
           }
           fprintf(ficgp,") t \"p%d%d\" ", k2,k);
           if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
           i=i+ncovmodel;
          }
        }
      }
      fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);
      }
      
     fclose(ficgp);
   }  /* end gnuplot */
   
   
   /*************** Moving average **************/
   void movingaverage(double agedeb, double fage,double agemin, double ***mobaverage){
   
     int i, cpt, cptcod;
       for (agedeb=agemin; agedeb<=fage; agedeb++)
         for (i=1; i<=nlstate;i++)
           for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
             mobaverage[(int)agedeb][i][cptcod]=0.;
      
       for (agedeb=agemin+4; agedeb<=fage; agedeb++){
         for (i=1; i<=nlstate;i++){
           for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
             for (cpt=0;cpt<=4;cpt++){
               mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod];
             }
             mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5;
           }
         }
       }
      
   }
   
   
   /************** Forecasting ******************/
   prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double agemin, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){
    
     int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
     int *popage;
     double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
     double *popeffectif,*popcount;
     double ***p3mat;
     char fileresf[FILENAMELENGTH];
   
    agelim=AGESUP;
   calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
    
     prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
   
    
     strcpy(fileresf,"f");
     strcat(fileresf,fileres);
     if((ficresf=fopen(fileresf,"w"))==NULL) {
       printf("Problem with forecast resultfile: %s\n", fileresf);
     }
     printf("Computing forecasting: result on file '%s' \n", fileresf);
   
     if (cptcoveff==0) ncodemax[cptcoveff]=1;
   
     if (mobilav==1) {
       mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
       movingaverage(agedeb, fage, agemin, mobaverage);
     }
   
     stepsize=(int) (stepm+YEARM-1)/YEARM;
     if (stepm<=12) stepsize=1;
    
     agelim=AGESUP;
    
     hstepm=1;
     hstepm=hstepm/stepm;
     yp1=modf(dateintmean,&yp);
     anprojmean=yp;
     yp2=modf((yp1*12),&yp);
     mprojmean=yp;
     yp1=modf((yp2*30.5),&yp);
     jprojmean=yp;
     if(jprojmean==0) jprojmean=1;
     if(mprojmean==0) jprojmean=1;
    
     fprintf(ficresf,"# Estimated date of observed prevalence: %.lf/%.lf/%.lf ",jprojmean,mprojmean,anprojmean);
    
     for(cptcov=1;cptcov<=i2;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
         k=k+1;
         fprintf(ficresf,"\n#******");
         for(j=1;j<=cptcoveff;j++) {
           fprintf(ficresf," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         }
         fprintf(ficresf,"******\n");
         fprintf(ficresf,"# StartingAge FinalAge");
         for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
        
        
         for (cpt=0; cpt<=(anproj2-anproj1);cpt++) {
           fprintf(ficresf,"\n");
           fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);  
          
           for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){
             nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
             nhstepm = nhstepm/hstepm;
            
             p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
             oldm=oldms;savm=savms;
             hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
          
             for (h=0; h<=nhstepm; h++){
               if (h==(int) (calagedate+YEARM*cpt)) {
                 fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);
               }
               for(j=1; j<=nlstate+ndeath;j++) {
                 kk1=0.;kk2=0;
                 for(i=1; i<=nlstate;i++) {              
                   if (mobilav==1)
                     kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];
                   else {
                     kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];
                   }
                  
                 }
                 if (h==(int)(calagedate+12*cpt)){
                   fprintf(ficresf," %.3f", kk1);
                          
                 }
               }
             }
             free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           }
         }
       }
     }
          
     if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
     fclose(ficresf);
   }
   /************** Forecasting ******************/
   populforecast(char fileres[], double anproj1,double mproj1,double jproj1,double agemin, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){
    
     int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
     int *popage;
     double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
     double *popeffectif,*popcount;
     double ***p3mat,***tabpop,***tabpopprev;
     char filerespop[FILENAMELENGTH];
   
   tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
    agelim=AGESUP;
   calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
    
     prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
   
    
     strcpy(filerespop,"pop");
     strcat(filerespop,fileres);
     if((ficrespop=fopen(filerespop,"w"))==NULL) {
       printf("Problem with forecast resultfile: %s\n", filerespop);
     }
     printf("Computing forecasting: result on file '%s' \n", filerespop);
   
     if (cptcoveff==0) ncodemax[cptcoveff]=1;
   
     if (mobilav==1) {
       mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
       movingaverage(agedeb, fage, agemin, mobaverage);
     }
   
     stepsize=(int) (stepm+YEARM-1)/YEARM;
     if (stepm<=12) stepsize=1;
    
     agelim=AGESUP;
    
     hstepm=1;
     hstepm=hstepm/stepm;
    
     if (popforecast==1) {
       if((ficpop=fopen(popfile,"r"))==NULL) {
         printf("Problem with population file : %s\n",popfile);exit(0);
       }
       popage=ivector(0,AGESUP);
       popeffectif=vector(0,AGESUP);
       popcount=vector(0,AGESUP);
      
       i=1;  
       while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF) i=i+1;
      
       imx=i;
       for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];
     }
   
     for(cptcov=1;cptcov<=i2;cptcov++){
      for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
         k=k+1;
         fprintf(ficrespop,"\n#******");
         for(j=1;j<=cptcoveff;j++) {
           fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         }
         fprintf(ficrespop,"******\n");
         fprintf(ficrespop,"# StartingAge FinalAge");
         for(j=1; j<=nlstate+ndeath;j++) fprintf(ficrespop," P.%d",j);
         if (popforecast==1)  fprintf(ficrespop," [Population]");
        
         for (cpt=0; cpt<=0;cpt++) {
           fprintf(ficrespop,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);  
          
           for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){
             nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
             nhstepm = nhstepm/hstepm;
            
             p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
             oldm=oldms;savm=savms;
             hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
          
             for (h=0; h<=nhstepm; h++){
               if (h==(int) (calagedate+YEARM*cpt)) {
                 fprintf(ficrespop,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);
               }
               for(j=1; j<=nlstate+ndeath;j++) {
                 kk1=0.;kk2=0;
                 for(i=1; i<=nlstate;i++) {              
                   if (mobilav==1)
                     kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];
                   else {
                     kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];
                   }
                 }
                 if (h==(int)(calagedate+12*cpt)){
                   tabpop[(int)(agedeb)][j][cptcod]=kk1;
                     /*fprintf(ficrespop," %.3f", kk1);
                       if (popforecast==1) fprintf(ficrespop," [%.f]", kk1*popeffectif[(int)agedeb+1]);*/
                 }
               }
               for(i=1; i<=nlstate;i++){
                 kk1=0.;
                   for(j=1; j<=nlstate;j++){
                     kk1= kk1+Tabpop[(int)(agedeb)][j][cptcod];
                   }
                     tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedate+12*cpt)*hstepm/YEARM*stepm-1)];
               }
   
   if (h==(int)(calagedate+12*cpt)) for(j=1; j<=nlstate;j++) fprintf(ficrespop," %.3f",tabpopprev[(int)(agedeb+1)][j][cptcod]);
             }
             free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           }
         }
    
     /******/
   
         for (cpt=1; cpt<=4;cpt++) {
           fprintf(ficrespop,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);  
           for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){
             nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
             nhstepm = nhstepm/hstepm;
            
             p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
             oldm=oldms;savm=savms;
             hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
             for (h=0; h<=nhstepm; h++){
               if (h==(int) (calagedate+YEARM*cpt)) {
                 fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);
               }
               for(j=1; j<=nlstate+ndeath;j++) {
                 kk1=0.;kk2=0;
                 for(i=1; i<=nlstate;i++) {              
                   kk1=kk1+p3mat[i][j][h]*tabpopprev[(int)agedeb+1][i][cptcod];    
                 }
                 if (h==(int)(calagedate+12*cpt)) fprintf(ficresf," %.3f", kk1);  
               }
             }
             free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           }
         }
      }
     }
    
     if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
     if (popforecast==1) {
       free_ivector(popage,0,AGESUP);
       free_vector(popeffectif,0,AGESUP);
       free_vector(popcount,0,AGESUP);
     }
     free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     fclose(ficrespop);
   }
   
   /***********************************************/
   /**************** Main Program *****************/
   /***********************************************/
   
   int main(int argc, char *argv[])
   {
   
     int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
     double agedeb, agefin,hf;
     double agemin=1.e20, agemax=-1.e20;
   
     double fret;
     double **xi,tmp,delta;
   
     double dum; /* Dummy variable */
     double ***p3mat;
     int *indx;
     char line[MAXLINE], linepar[MAXLINE];
     char title[MAXLINE];
     char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];
     char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
    
     char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];
   
     char filerest[FILENAMELENGTH];
     char fileregp[FILENAMELENGTH];
     char popfile[FILENAMELENGTH];
     char path[80],pathc[80],pathcd[80],pathtot[80],model[20];
     int firstobs=1, lastobs=10;
     int sdeb, sfin; /* Status at beginning and end */
     int c,  h , cpt,l;
     int ju,jl, mi;
     int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
     int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
     int mobilav=0,popforecast=0;
     int hstepm, nhstepm;
     double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2;
   
     double bage, fage, age, agelim, agebase;
     double ftolpl=FTOL;
     double **prlim;
     double *severity;
     double ***param; /* Matrix of parameters */
     double  *p;
     double **matcov; /* Matrix of covariance */
     double ***delti3; /* Scale */
     double *delti; /* Scale */
     double ***eij, ***vareij;
     double **varpl; /* Variances of prevalence limits by age */
     double *epj, vepp;
     double kk1, kk2;
     double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
    
   
     char version[80]="Imach version 0.7, February 2002, INED-EUROREVES ";
     char *alph[]={"a","a","b","c","d","e"}, str[4];
   
   
     char z[1]="c", occ;
   #include <sys/time.h>
   #include <time.h>
     char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];
    
     /* long total_usecs;
     struct timeval start_time, end_time;
    
     gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
   
   
     printf("\n%s",version);
     if(argc <=1){
       printf("\nEnter the parameter file name: ");
       scanf("%s",pathtot);
     }
     else{
       strcpy(pathtot,argv[1]);
     }
     /*if(getcwd(pathcd, 80)!= NULL)printf ("Error pathcd\n");*/
     /*cygwin_split_path(pathtot,path,optionfile);
       printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
     /* cutv(path,optionfile,pathtot,'\\');*/
   
     split(pathtot,path,optionfile,optionfilext,optionfilefiname);
      printf("pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
     chdir(path);
     replace(pathc,path);
   
   /*-------- arguments in the command line --------*/
   
     strcpy(fileres,"r");
     strcat(fileres, optionfilefiname);
     strcat(fileres,".txt");    /* Other files have txt extension */
   
     /*---------arguments file --------*/
   
     if((ficpar=fopen(optionfile,"r"))==NULL)    {
       printf("Problem with optionfile %s\n",optionfile);
       goto end;
     }
   
     strcpy(filereso,"o");
     strcat(filereso,fileres);
     if((ficparo=fopen(filereso,"w"))==NULL) {
       printf("Problem with Output resultfile: %s\n", filereso);goto end;
     }
   
     /* Reads comments: lines beginning with '#' */
     while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
   
     fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
     printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt,model);
     fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt,model);
   while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
    
      
     covar=matrix(0,NCOVMAX,1,n);
     cptcovn=0;
     if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
   
     ncovmodel=2+cptcovn;
     nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
    
     /* Read guess parameters */
     /* Reads comments: lines beginning with '#' */
     while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
    
     param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
       for(i=1; i <=nlstate; i++)
       for(j=1; j <=nlstate+ndeath-1; j++){
         fscanf(ficpar,"%1d%1d",&i1,&j1);
         fprintf(ficparo,"%1d%1d",i1,j1);
         printf("%1d%1d",i,j);
         for(k=1; k<=ncovmodel;k++){
           fscanf(ficpar," %lf",&param[i][j][k]);
           printf(" %lf",param[i][j][k]);
           fprintf(ficparo," %lf",param[i][j][k]);
         }
         fscanf(ficpar,"\n");
         printf("\n");
         fprintf(ficparo,"\n");
       }
    
       npar= (nlstate+ndeath-1)*nlstate*ncovmodel;
   
     p=param[1][1];
    
     /* Reads comments: lines beginning with '#' */
     while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
   
     delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
     delti=vector(1,npar); /* Scale of each paramater (output from hesscov) */
     for(i=1; i <=nlstate; i++){
       for(j=1; j <=nlstate+ndeath-1; j++){
       fscanf(ficpar,"%1d%1d",&i1,&j1);        fscanf(ficpar,"%1d%1d",&i1,&j1);
       printf("%1d%1d",i,j);        printf("%1d%1d",i,j);
       fprintf(ficparo,"%1d%1d",i1,j1);        fprintf(ficparo,"%1d%1d",i1,j1);
Line 1862  split(pathtot, path,optionfile); Line 2652  split(pathtot, path,optionfile);
     tab=ivector(1,NCOVMAX);      tab=ivector(1,NCOVMAX);
     ncodemax=ivector(1,8);      ncodemax=ivector(1,8);
   
     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)) {
                 
Line 1884  split(pathtot, path,optionfile); Line 2674  split(pathtot, path,optionfile);
           cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);            cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
         }          }
         num[i]=atol(stra);          num[i]=atol(stra);
          
         /*printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/          /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
             printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
   
         i=i+1;          i=i+1;
       }        }
     }      }
       /* printf("ii=%d", ij);
     /*scanf("%d",i);*/         scanf("%d",i);*/
   imx=i-1; /* Number of individuals */    imx=i-1; /* Number of individuals */
   
     /* for (i=1; i<=imx; i++){
       if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;
       if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
       if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
       }
   
       for (i=1; i<=imx; i++)
       if (covar[1][i]==0) printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/
   
   /* Calculation of the number of parameter from char model*/    /* Calculation of the number of parameter from char model*/
   Tvar=ivector(1,15);    Tvar=ivector(1,15);
   Tprod=ivector(1,15);    Tprod=ivector(1,15);
Line 1975  split(pathtot, path,optionfile); Line 2775  split(pathtot, path,optionfile);
     }      }
     /*-calculation of age at interview from date of interview and age at death -*/      /*-calculation of age at interview from date of interview and age at death -*/
     agev=matrix(1,maxwav,1,imx);      agev=matrix(1,maxwav,1,imx);
   
      for (i=1; i<=imx; i++)
        for(m=2; (m<= maxwav); m++)
          if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
            anint[m][i]=9999;
            s[m][i]=-1;
          }
         
     for (i=1; i<=imx; i++)  {      for (i=1; i<=imx; i++)  {
       agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);        agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
Line 2031  printf("Total number of individuals= %d, Line 2838  printf("Total number of individuals= %d,
     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);
     free_vector(annais,1,n);      free_vector(annais,1,n);
     free_matrix(mint,1,maxwav,1,n);      /* free_matrix(mint,1,maxwav,1,n);
     free_matrix(anint,1,maxwav,1,n);         free_matrix(anint,1,maxwav,1,n);*/
     free_vector(moisdc,1,n);      free_vector(moisdc,1,n);
     free_vector(andc,1,n);      free_vector(andc,1,n);
   
Line 2076  printf("Total number of individuals= %d, Line 2883  printf("Total number of individuals= %d,
         
    /* Calculates basic frequencies. Computes observed prevalence at single age     /* Calculates basic frequencies. Computes observed prevalence at single age
        and prints on file fileres'p'. */         and prints on file fileres'p'. */
   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax);  
   
     pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */  
     oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */  
     newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */  
     savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */  
     oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */  
      
     /* For Powell, parameters are in a vector p[] starting at p[1]  
        so we point p on param[1][1] so that p[1] maps on param[1][1][1] */  
     p=param[1][1]; /* *(*(*(param +1)+1)+0) */  
   
     if(mle==1){  
     mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);  
     }  
      
     /*--------- results files --------------*/  
     fprintf(ficres,"\ntitle=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, mle,weightopt,model);  
      
    jk=1;  
    fprintf(ficres,"# Parameters\n");  
    printf("# Parameters\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(ficres,"%1d%1d ",i,k);  
            for(j=1; j <=ncovmodel; j++){  
              printf("%f ",p[jk]);  
              fprintf(ficres,"%f ",p[jk]);  
              jk++;  
            }  
            printf("\n");  
            fprintf(ficres,"\n");  
          }  
      }  
    }  
  if(mle==1){  
     /* Computing hessian and covariance matrix */  
     ftolhess=ftol; /* Usually correct */  
     hesscov(matcov, p, npar, delti, ftolhess, func);  
  }  
     fprintf(ficres,"# Scales\n");  
     printf("# Scales\n");  
      for(i=1,jk=1; i <=nlstate; i++){  
       for(j=1; j <=nlstate+ndeath; j++){  
         if (j!=i) {  
           fprintf(ficres,"%1d%1d",i,j);  
           printf("%1d%1d",i,j);  
           for(k=1; k<=ncovmodel;k++){  
             printf(" %.5e",delti[jk]);  
             fprintf(ficres," %.5e",delti[jk]);  
             jk++;  
           }  
           printf("\n");  
           fprintf(ficres,"\n");  
         }  
       }  
       }  
      
     k=1;  
     fprintf(ficres,"# Covariance\n");  
     printf("# Covariance\n");  
     for(i=1;i<=npar;i++){  
       /*  if (k>nlstate) k=1;  
       i1=(i-1)/(ncovmodel*nlstate)+1;  
       fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);  
       printf("%s%d%d",alph[k],i1,tab[i]);*/  
       fprintf(ficres,"%3d",i);  
       printf("%3d",i);  
       for(j=1; j<=i;j++){  
         fprintf(ficres," %.5e",matcov[i][j]);  
         printf(" %.5e",matcov[i][j]);  
       }  
       fprintf(ficres,"\n");  
       printf("\n");  
       k++;  
     }  
      
     while((c=getc(ficpar))=='#' && c!= EOF){  
       ungetc(c,ficpar);  
       fgets(line, MAXLINE, ficpar);  
       puts(line);  
       fputs(line,ficparo);  
     }  
     ungetc(c,ficpar);  
    
     fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemax, &bage, &fage);  
      
     if (fage <= 2) {  
       bage = agemin;  
       fage = agemax;  
     }  
   
     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\n",agemin,agemax,bage,fage);  
   
      
 /*------------ gnuplot -------------*/  
 chdir(pathcd);  
   if((ficgp=fopen("graph.plt","w"))==NULL) {  
     printf("Problem with file graph.gp");goto end;  
   }  
 #ifdef windows  
   fprintf(ficgp,"cd \"%s\" \n",pathc);  
 #endif  
 m=pow(2,cptcoveff);  
    
  /* 1eme*/  
   for (cpt=1; cpt<= nlstate ; cpt ++) {  
    for (k1=1; k1<= m ; k1 ++) {  
   
 #ifdef windows  
     fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",agemin,fage,fileres,k1-1,k1-1);  
 #endif  
 #ifdef unix  
 fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",agemin,fage,fileres);  
 #endif  
   
 for (i=1; i<= nlstate ; i ++) {  
   if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");  
   else fprintf(ficgp," \%%*lf (\%%*lf)");  
 }  
     fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1);  
     for (i=1; i<= nlstate ; i ++) {  
   if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");  
   else fprintf(ficgp," \%%*lf (\%%*lf)");  
 }  
   fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"vpl%s\" every :::%d::%d u 1:($2-2*$3) \"\%%lf",fileres,k1-1,k1-1);  
      for (i=1; i<= nlstate ; i ++) {  
   if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");  
   else fprintf(ficgp," \%%*lf (\%%*lf)");  
 }    
      fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));  
 #ifdef unix  
 fprintf(ficgp,"\nset ter gif small size 400,300");  
 #endif  
 fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);  
    }  
   }  
   /*2 eme*/  
   
   for (k1=1; k1<= m ; k1 ++) {  
     fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",agemin,fage);  
      
     for (i=1; i<= nlstate+1 ; i ++) {  
       k=2*i;  
       fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:2 \"\%%lf",fileres,k1-1,k1-1);  
       for (j=1; j<= nlstate+1 ; j ++) {  
   if (j==i) fprintf(ficgp," \%%lf (\%%lf)");  
   else fprintf(ficgp," \%%*lf (\%%*lf)");  
 }    
       if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");  
       else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);  
     fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",fileres,k1-1,k1-1);  
       for (j=1; j<= nlstate+1 ; j ++) {  
         if (j==i) fprintf(ficgp," \%%lf (\%%lf)");  
         else fprintf(ficgp," \%%*lf (\%%*lf)");  
 }    
       fprintf(ficgp,"\" t\"\" w l 0,");  
      fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",fileres,k1-1,k1-1);  
       for (j=1; j<= nlstate+1 ; j ++) {  
   if (j==i) fprintf(ficgp," \%%lf (\%%lf)");  
   else fprintf(ficgp," \%%*lf (\%%*lf)");  
 }    
       if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");  
       else fprintf(ficgp,"\" t\"\" w l 0,");  
     }  
     fprintf(ficgp,"\nset out \"e%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),k1);  
   }  
    
   /*3eme*/  
   
   for (k1=1; k1<= m ; k1 ++) {  
     for (cpt=1; cpt<= nlstate ; cpt ++) {  
       k=2+nlstate*(cpt-1);  
       fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",agemin,fage,fileres,k1-1,k1-1,k,cpt);  
       for (i=1; i< nlstate ; i ++) {  
         fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+i,cpt,i+1);  
       }  
       fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);  
     }  
   }  
    
   /* CV preval stat */  
   for (k1=1; k1<= m ; k1 ++) {  
     for (cpt=1; cpt<nlstate ; cpt ++) {  
       k=3;  
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",agemin,agemax,fileres,k1,k+cpt+1,k+1);  
       for (i=1; i< nlstate ; i ++)  
         fprintf(ficgp,"+$%d",k+i+1);  
       fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);  
        
       l=3+(nlstate+ndeath)*cpt;  
       fprintf(ficgp,",\"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",fileres,k1,l+cpt+1,l+1);  
       for (i=1; i< nlstate ; i ++) {  
         l=3+(nlstate+ndeath)*cpt;  
         fprintf(ficgp,"+$%d",l+i+1);  
       }  
       fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);    
       fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);  
     }  
   }  
   
   /* proba elementaires */     
      
       pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
       oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
       newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
       savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
       oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
        
       /* For Powell, parameters are in a vector p[] starting at p[1]
          so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
       p=param[1][1]; /* *(*(*(param +1)+1)+0) */
   
       if(mle==1){
       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 ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, mle,weightopt,model);
    
   
      jk=1;
      fprintf(ficres,"# Parameters\n");
      printf("# Parameters\n");
    for(i=1,jk=1; i <=nlstate; i++){     for(i=1,jk=1; i <=nlstate; i++){
     for(k=1; k <=(nlstate+ndeath); k++){       for(k=1; k <=(nlstate+ndeath); k++){
       if (k != i) {         if (k != i)
         for(j=1; j <=ncovmodel; j++){           {
           /*fprintf(ficgp,"%s%1d%1d=%f ",alph[j],i,k,p[jk]);*/             printf("%d%d ",i,k);
           /*fprintf(ficgp,"%s",alph[1]);*/             fprintf(ficres,"%1d%1d ",i,k);
           fprintf(ficgp,"p%d=%f ",jk,p[jk]);             for(j=1; j <=ncovmodel; j++){
           jk++;               printf("%f ",p[jk]);
           fprintf(ficgp,"\n");               fprintf(ficres,"%f ",p[jk]);
                jk++;
              }
              printf("\n");
              fprintf(ficres,"\n");
            }
        }
      }
    if(mle==1){
       /* Computing hessian and covariance matrix */
       ftolhess=ftol; /* Usually correct */
       hesscov(matcov, p, npar, delti, ftolhess, func);
    }
       fprintf(ficres,"# Scales\n");
       printf("# Scales\n");
        for(i=1,jk=1; i <=nlstate; i++){
         for(j=1; j <=nlstate+ndeath; j++){
           if (j!=i) {
             fprintf(ficres,"%1d%1d",i,j);
             printf("%1d%1d",i,j);
             for(k=1; k<=ncovmodel;k++){
               printf(" %.5e",delti[jk]);
               fprintf(ficres," %.5e",delti[jk]);
               jk++;
             }
             printf("\n");
             fprintf(ficres,"\n");
         }          }
       }        }
        }
      
       k=1;
       fprintf(ficres,"# Covariance\n");
       printf("# Covariance\n");
       for(i=1;i<=npar;i++){
         /*  if (k>nlstate) k=1;
         i1=(i-1)/(ncovmodel*nlstate)+1;
         fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);
         printf("%s%d%d",alph[k],i1,tab[i]);*/
         fprintf(ficres,"%3d",i);
         printf("%3d",i);
         for(j=1; j<=i;j++){
           fprintf(ficres," %.5e",matcov[i][j]);
           printf(" %.5e",matcov[i][j]);
         }
         fprintf(ficres,"\n");
         printf("\n");
         k++;
       }
      
       while((c=getc(ficpar))=='#' && c!= EOF){
         ungetc(c,ficpar);
         fgets(line, MAXLINE, ficpar);
         puts(line);
         fputs(line,ficparo);
     }      }
       ungetc(c,ficpar);
    
       fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemax, &bage, &fage);
      
       if (fage <= 2) {
         bage = agemin;
         fage = agemax;
     }      }
      
       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\n",agemin,agemax,bage,fage);
       fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
    
       while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
    
     fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2);
     fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
    fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
        
     while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
    
   
   for(jk=1; jk <=m; jk++) {     dateprev1=anprev1+mprev1/12.+jprev1/365.;
   fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",agemin,agemax);     dateprev2=anprev2+mprev2/12.+jprev2/365.;
    i=1;  
    for(k2=1; k2<=nlstate; k2++) {    fscanf(ficpar,"pop_based=%d\n",&popbased);
      k3=i;     fprintf(ficparo,"pop_based=%d\n",popbased);  
      for(k=1; k<=(nlstate+ndeath); k++) {     fprintf(ficres,"pop_based=%d\n",popbased);  
        if (k != k2){  
         fprintf(ficgp," exp(p%d+p%d*x",i,i+1);    while((c=getc(ficpar))=='#' && c!= EOF){
 ij=1;      ungetc(c,ficpar);
         for(j=3; j <=ncovmodel; j++) {      fgets(line, MAXLINE, ficpar);
           if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {      puts(line);
             fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);      fputs(line,ficparo);
             ij++;  
           }  
           else  
           fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);  
         }  
           fprintf(ficgp,")/(1");  
          
         for(k1=1; k1 <=nlstate; k1++){    
           fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);  
 ij=1;  
           for(j=3; j <=ncovmodel; j++){  
           if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {  
             fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);  
             ij++;  
           }  
           else  
             fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);  
           }  
           fprintf(ficgp,")");  
         }  
         fprintf(ficgp,") t \"p%d%d\" ", k2,k);  
         if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");  
         i=i+ncovmodel;  
        }  
      }  
    }  
    fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);  
   }    }
        ungetc(c,ficpar);
   fclose(ficgp);    fscanf(ficpar,"popforecast=%d popfile=%s starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mov_average=%d\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilav);
      fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);
 chdir(path);  fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);
     free_matrix(agev,1,maxwav,1,imx);  
     free_ivector(wav,1,imx);  
     free_imatrix(dh,1,lastpass-firstpass+1,1,imx);  
     free_imatrix(mw,1,lastpass-firstpass+1,1,imx);  
      
     free_imatrix(s,1,maxwav+1,1,n);  
      
      
     free_ivector(num,1,n);  
     free_vector(agedc,1,n);  
     free_vector(weight,1,n);  
     /*free_matrix(covar,1,NCOVMAX,1,n);*/  
     fclose(ficparo);  
     fclose(ficres);  
     /*  }*/  
      
    /*________fin mle=1_________*/  
      
   
     freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
     /* No more information from the sample is required now */  
   /*------------ gnuplot -------------*/
    printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, agemin,agemax,fage, pathc,p);
    
   /*------------ free_vector  -------------*/
    chdir(path);
    
    free_ivector(wav,1,imx);
    free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
    free_imatrix(mw,1,lastpass-firstpass+1,1,imx);  
    free_ivector(num,1,n);
    free_vector(agedc,1,n);
    /*free_matrix(covar,1,NCOVMAX,1,n);*/
    fclose(ficparo);
    fclose(ficres);
    
   /* 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);      ungetc(c,ficpar);
Line 2375  chdir(path); Line 3054  chdir(path);
   fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);    fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
 /*--------- index.htm --------*/  /*--------- index.htm --------*/
   
   strcpy(optionfilehtm,optionfile);    printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm);
   strcat(optionfilehtm,".htm");  
   if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {  
     printf("Problem with %s \n",optionfilehtm);goto end;  
   }  
   
  fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.64b </font> <hr size=\"2\" color=\"#EC5E5E\">  
 Titre=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>  
 Total number of observations=%d <br>  
 Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>  
 <hr  size=\"2\" color=\"#EC5E5E\">  
 <li>Outputs files<br><br>\n  
         - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n  
 - Estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>  
         - Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>  
         - Transition probabilities: <a href=\"pij%s\">pij%s</a><br>  
         - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>  
         - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>  
         - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>  
         - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>  
         - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br><br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);  
   
  fprintf(fichtm," <li>Graphs</li><p>");  
   
  m=cptcoveff;  
  if (cptcovn < 1) {m=1;ncodemax[1]=1;}  
   
  j1=0;  
  for(k1=1; k1<=m;k1++){  
    for(i1=1; i1<=ncodemax[k1];i1++){  
        j1++;  
        if (cptcovn > 0) {  
          fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");  
          for (cpt=1; cpt<=cptcoveff;cpt++)  
            fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[j1][cpt]]);  
          fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");  
        }  
        fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>  
 <img src=\"pe%s%d.gif\">",strtok(optionfile, "."),j1,strtok(optionfile, "."),j1);      
        for(cpt=1; cpt<nlstate;cpt++){  
          fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.gif<br>  
 <img src=\"p%s%d%d.gif\">",strtok(optionfile, "."),cpt,j1,strtok(optionfile, "."),cpt,j1);  
        }  
     for(cpt=1; cpt<=nlstate;cpt++) {  
        fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident  
 interval) in state (%d): v%s%d%d.gif <br>  
 <img src=\"v%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,j1,strtok(optionfile, "."),cpt,j1);    
      }  
      for(cpt=1; cpt<=nlstate;cpt++) {  
         fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.gif <br>  
 <img src=\"exp%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,j1,strtok(optionfile, "."),cpt,j1);  
      }  
      fprintf(fichtm,"\n<br>- Total life expectancy by age and  
 health expectancies in states (1) and (2): e%s%d.gif<br>  
 <img src=\"e%s%d.gif\">",strtok(optionfile, "."),j1,strtok(optionfile, "."),j1);  
 fprintf(fichtm,"\n</body>");  
    }  
  }  
 fclose(fichtm);  
   
    
   /*--------------- Prevalence limit --------------*/    /*--------------- Prevalence limit --------------*/
     
   strcpy(filerespl,"pl");    strcpy(filerespl,"pl");
Line 2480  fclose(fichtm); Line 3102  fclose(fichtm);
       }        }
     }      }
   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);
Line 2489  fclose(fichtm); Line 3112  fclose(fichtm);
   printf("Computing pij: result on file '%s' \n", filerespij);    printf("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 */
Line 2528  fclose(fichtm); Line 3151  fclose(fichtm);
     }      }
   }    }
   
     /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/
   
   fclose(ficrespij);    fclose(ficrespij);
   
   
     /*---------- Forecasting ------------------*/
     if(stepm == 1) {
       prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);
   populforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);
       free_matrix(mint,1,maxwav,1,n);
       free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);
       free_vector(weight,1,n);}
     else{
       erreur=108;
       printf("Error %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d\n", erreur, stepm);
     }
    
   
   /*---------- Health expectancies and variances ------------*/    /*---------- Health expectancies and variances ------------*/
   
   strcpy(filerest,"t");    strcpy(filerest,"t");
Line 2578  fclose(fichtm); Line 3217  fclose(fichtm);
       evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k);          evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k);  
       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(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);         varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);
           
   
    
       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");
          
       hf=1;        hf=1;
       if (stepm >= YEARM) hf=stepm/YEARM;        if (stepm >= YEARM) hf=stepm/YEARM;
       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) {
             for(i=1; i<=nlstate;i++)
               prlim[i][i]=probs[(int)age][i][k];
           }
          
         fprintf(ficrest," %.0f",age);          fprintf(ficrest," %.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++) {
Line 2607  fclose(fichtm); Line 3253  fclose(fichtm);
       }        }
     }      }
   }    }
          
  fclose(ficreseij);    fclose(ficreseij);
  fclose(ficresvij);    fclose(ficresvij);
   fclose(ficrest);    fclose(ficrest);
   fclose(ficpar);    fclose(ficpar);
   free_vector(epj,1,nlstate+1);    free_vector(epj,1,nlstate+1);
   /*  scanf("%d ",i); */   
   
   /*------- Variance limit prevalence------*/      /*------- Variance limit prevalence------*/  
   
 strcpy(fileresvpl,"vpl");    strcpy(fileresvpl,"vpl");
   strcat(fileresvpl,fileres);    strcat(fileresvpl,fileres);
   if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {    if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
     printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);      printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);
Line 2625  strcpy(fileresvpl,"vpl"); Line 3270  strcpy(fileresvpl,"vpl");
   }    }
   printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);    printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);
   
  k=0;    k=0;
  for(cptcov=1;cptcov<=i1;cptcov++){    for(cptcov=1;cptcov<=i1;cptcov++){
    for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){      for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
      k=k+1;        k=k+1;
      fprintf(ficresvpl,"\n#****** ");        fprintf(ficresvpl,"\n#****** ");
      for(j=1;j<=cptcoveff;j++)        for(j=1;j<=cptcoveff;j++)
        fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);          fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
      fprintf(ficresvpl,"******\n");        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);
    }      }
  }   }
   
   fclose(ficresvpl);    fclose(ficresvpl);
Line 2653  strcpy(fileresvpl,"vpl"); Line 3298  strcpy(fileresvpl,"vpl");
   free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);    free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);    free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);    free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
     
   free_matrix(matcov,1,npar,1,npar);    free_matrix(matcov,1,npar,1,npar);
   free_vector(delti,1,npar);    free_vector(delti,1,npar);
      free_matrix(agev,1,maxwav,1,imx);
   free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);    free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
   
   printf("End of Imach\n");    if(erreur >0)
       printf("End of Imach with error %d\n",erreur);
     else   printf("End of Imach\n");
   /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */    /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */
     
   /* printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);*/    /* printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);*/
   /*printf("Total time was %d uSec.\n", total_usecs);*/    /*printf("Total time was %d uSec.\n", total_usecs);*/
   /*------ End -----------*/    /*------ End -----------*/
   
   
  end:   end:
 #ifdef windows  #ifdef windows
  chdir(pathcd);    /* chdir(pathcd);*/
 #endif  #endif
  /*system("wgnuplot graph.plt");*/   /*system("wgnuplot graph.plt");*/
  system("../gp37mgw/wgnuplot graph.plt");   /*system("../gp37mgw/wgnuplot graph.plt");*/
    /*system("cd ../gp37mgw");*/
    /* system("..\\gp37mgw\\wgnuplot graph.plt");*/
    strcpy(plotcmd,GNUPLOTPROGRAM);
    strcat(plotcmd," ");
    strcat(plotcmd,optionfilegnuplot);
    system(plotcmd);
   
 #ifdef windows  #ifdef windows
   while (z[0] != 'q') {    while (z[0] != 'q') {
     chdir(pathcd);      chdir(path);
     printf("\nType e to edit output files, c to start again, and q for exiting: ");      printf("\nType e to edit output files, c to start again, and q for exiting: ");
     scanf("%s",z);      scanf("%s",z);
     if (z[0] == 'c') system("./imach");      if (z[0] == 'c') system("./imach");
Removed from v.1.10  
changed lines
  Added in v.1.27

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