Diff for /imach/src/imach.c between versions 1.1.1.1 and 1.42

version 1.1.1.1, 2000/12/28 18:49:56 version 1.42, 2002/05/21 18:44:41
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 Likekilhood 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     probability 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
   Also this programme outputs the covariance matrix of the parameters but also    conditional to the observed state i at age x. The delay 'h' can be
   of the life expectancies. It also computes the prevalence limits.     split into an exact number (nh*stepm) of unobserved intermediate
       states. This elementary transition (by month or quarter trimester,
   Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).    semester or year) is model as a multinomial logistic.  The hPx
            Institut national d'études démographiques, Paris.    matrix is simply the matrix product of nh*stepm elementary matrices
   This software have been partly granted by Euro-REVES, a concerted action    and the contribution of each individual to the likelihood is simply
   from the European Union.    hPijx.
   It is copyrighted identically to a GNU software product, ie programme and  
   software can be distributed freely for non commercial use. Latest version    Also this programme outputs the covariance matrix of the parameters but also
   can be accessed at http://euroreves.ined.fr/imach .    of the life expectancies. It also computes the prevalence limits.
   **********************************************************************/   
      Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).
 #include <math.h>             Institut national d'études démographiques, Paris.
 #include <stdio.h>    This software have been partly granted by Euro-REVES, a concerted action
 #include <stdlib.h>    from the European Union.
 #include <unistd.h>    It is copyrighted identically to a GNU software product, ie programme and
     software can be distributed freely for non commercial use. Latest version
 #define MAXLINE 256    can be accessed at http://euroreves.ined.fr/imach .
 #define FILENAMELENGTH 80    **********************************************************************/
 /*#define DEBUG*/   
 /*#define win*/  #include <math.h>
   #include <stdio.h>
 #define MAXPARM 30 /* Maximum number of parameters for the optimization */  #include <stdlib.h>
 #define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncov */  #include <unistd.h>
   
 #define NINTERVMAX 8  #define MAXLINE 256
 #define NLSTATEMAX 8 /* Maximum number of live states (for func) */  #define GNUPLOTPROGRAM "gnuplot"
 #define NDEATHMAX 8 /* Maximum number of dead states (for func) */  /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
 #define NCOVMAX 8 /* Maximum number of covariates */  #define FILENAMELENGTH 80
 #define MAXN 20000  /*#define DEBUG*/
 #define YEARM 12. /* Number of months per year */  #define windows
 #define AGESUP 130  #define GLOCK_ERROR_NOPATH              -1      /* empty path */
 #define AGEBASE 40  #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */
   
   #define MAXPARM 30 /* Maximum number of parameters for the optimization */
 int nvar;  #define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */
   
 int npar=NPARMAX;  #define NINTERVMAX 8
 int nlstate=2; /* Number of live states */  #define NLSTATEMAX 8 /* Maximum number of live states (for func) */
 int ndeath=1; /* Number of dead states */  #define NDEATHMAX 8 /* Maximum number of dead states (for func) */
 int ncov;     /* Total number of covariables including constant a12*1 +b12*x ncov=2 */  #define NCOVMAX 8 /* Maximum number of covariates */
   #define MAXN 20000
 int *wav; /* Number of waves for this individuual 0 is possible */  #define YEARM 12. /* Number of months per year */
 int maxwav; /* Maxim number of waves */  #define AGESUP 130
 int mle, weightopt;  #define AGEBASE 40
 int **mw; /* mw[mi][i] is number of the mi wave for this individual */  
 int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */  
 double **oldm, **newm, **savm; /* Working pointers to matrices */  int erreur; /* Error number */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */  int nvar;
 FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest;  int cptcovn, cptcovage=0, cptcoveff=0,cptcov;
 FILE *ficgp, *fichtm;  int npar=NPARMAX;
   int nlstate=2; /* Number of live states */
 #define NR_END 1  int ndeath=1; /* Number of dead states */
 #define FREE_ARG char*  int ncovmodel, ncovcol;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
 #define FTOL 1.0e-10  int popbased=0;
   
 #define NRANSI   int *wav; /* Number of waves for this individuual 0 is possible */
 #define ITMAX 200   int maxwav; /* Maxim number of waves */
   int jmin, jmax; /* min, max spacing between 2 waves */
 #define TOL 2.0e-4   int mle, weightopt;
   int **mw; /* mw[mi][i] is number of the mi wave for this individual */
 #define CGOLD 0.3819660   int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */
 #define ZEPS 1.0e-10   double jmean; /* Mean space between 2 waves */
 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);   double **oldm, **newm, **savm; /* Working pointers to matrices */
   double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
 #define GOLD 1.618034   FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
 #define GLIMIT 100.0   FILE *ficgp,*ficresprob,*ficpop;
 #define TINY 1.0e-20   FILE *ficreseij;
     char filerese[FILENAMELENGTH];
 static double maxarg1,maxarg2;   FILE  *ficresvij;
 #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)>(maxarg2)? (maxarg1):(maxarg2))    char fileresv[FILENAMELENGTH];
 #define FMIN(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)<(maxarg2)? (maxarg1):(maxarg2))   FILE  *ficresvpl;
      char fileresvpl[FILENAMELENGTH];
 #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a))  
 #define rint(a) floor(a+0.5)  #define NR_END 1
   #define FREE_ARG char*
 static double sqrarg;  #define FTOL 1.0e-10
 #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 :sqrarg*sqrarg)  
 #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}   #define NRANSI
   #define ITMAX 200
 int imx;   
 int stepm;  #define TOL 2.0e-4
 /* Stepm, step in month: minimum step interpolation*/  
   #define CGOLD 0.3819660
 int m,nb;  #define ZEPS 1.0e-10
 int *num, firstpass=0, lastpass=2,*cod;  #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
 double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;  
 double **pmmij;  #define GOLD 1.618034
   #define GLIMIT 100.0
 double *weight;  #define TINY 1.0e-20
 int **s; /* Status */  
 double *agedc, **covar, idx;  static double 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))
 double ftol=FTOL; /* Tolerance for computing Max Likelihood */   
 double ftolhess; /* Tolerance for computing hessian */  #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a))
   #define rint(a) floor(a+0.5)
   
 /******************************************/  static double sqrarg;
   #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 :sqrarg*sqrarg)
 void replace(char *s, char*t)  #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
 {  
   int i;  int imx;
   int lg=20;  int stepm;
   i=0;  /* Stepm, step in month: minimum step interpolation*/
   lg=strlen(t);  
   for(i=0; i<= lg; i++) {  int estepm;
     (s[i] = t[i]);  /* Estepm, step in month to interpolate survival function in order to approximate Life Expectancy*/
     if (t[i]== '\\') s[i]='/';  
   }  int m,nb;
 }  int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;
 void cut(char *u,char *v, char*t)  double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
 {  double **pmmij, ***probs, ***mobaverage;
   int i,lg,j,p;  double dateintmean=0;
   i=0;  
   for(j=0; j<=strlen(t); j++) {  double *weight;
     if(t[j]=='\\') p=j;  int **s; /* Status */
   }  double *agedc, **covar, idx;
   int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff;
   lg=strlen(t);  
   for(j=0; j<p; j++) {  double ftol=FTOL; /* Tolerance for computing Max Likelihood */
     (u[j] = t[j]);  double ftolhess; /* Tolerance for computing hessian */
     u[p]='\0';  
   }  /**************** split *************************/
   static  int split( char *path, char *dirc, char *name, char *ext, char *finame )
   for(j=0; j<= lg; j++) {  {
     if (j>=(p+1))(v[j-p-1] = t[j]);     char *s;                             /* pointer */
   }     int  l1, l2;                         /* length counters */
 }  
      l1 = strlen( path );                 /* length of path */
 /********************** nrerror ********************/     if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
   #ifdef windows
 void nrerror(char error_text[])     s = strrchr( path, '\\' );           /* find last / */
 {  #else
   fprintf(stderr,"ERREUR ...\n");     s = strrchr( path, '/' );            /* find last / */
   fprintf(stderr,"%s\n",error_text);  #endif
   exit(1);     if ( s == NULL ) {                   /* no directory, so use current */
 }  #if     defined(__bsd__)                /* get current working directory */
 /*********************** vector *******************/        extern char       *getwd( );
 double *vector(int nl, int nh)  
 {        if ( getwd( dirc ) == NULL ) {
   double *v;  #else
   v=(double *) malloc((size_t)((nh-nl+1+NR_END)*sizeof(double)));        extern char       *getcwd( );
   if (!v) nrerror("allocation failure in vector");  
   return v-nl+NR_END;        if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {
 }  #endif
            return( GLOCK_ERROR_GETCWD );
 /************************ free vector ******************/        }
 void free_vector(double*v, int nl, int nh)        strcpy( name, path );             /* we've got it */
 {     } else {                             /* strip direcotry from path */
   free((FREE_ARG)(v+nl-NR_END));        s++;                              /* after this, the filename */
 }        l2 = strlen( s );                 /* length of filename */
         if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH );
 /************************ivector *******************************/        strcpy( name, s );                /* save file name */
 int *ivector(long nl,long nh)        strncpy( dirc, path, l1 - l2 );   /* now the directory */
 {        dirc[l1-l2] = 0;                  /* add zero */
   int *v;     }
   v=(int *) malloc((size_t)((nh-nl+1+NR_END)*sizeof(int)));     l1 = strlen( dirc );                 /* length of directory */
   if (!v) nrerror("allocation failure in ivector");  #ifdef windows
   return v-nl+NR_END;     if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
 }  #else
      if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }
 /******************free ivector **************************/  #endif
 void free_ivector(int *v, long nl, long nh)     s = strrchr( name, '.' );            /* find last / */
 {     s++;
   free((FREE_ARG)(v+nl-NR_END));     strcpy(ext,s);                       /* save extension */
 }     l1= strlen( name);
      l2= strlen( s)+1;
 /******************* imatrix *******************************/     strncpy( finame, name, l1-l2);
 int **imatrix(long nrl, long nrh, long ncl, long nch)      finame[l1-l2]= 0;
      /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */      return( 0 );                         /* we're done */
 {   }
   long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;   
   int **m;   
     /******************************************/
   /* allocate pointers to rows */   
   m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));   void replace(char *s, char*t)
   if (!m) nrerror("allocation failure 1 in matrix()");   {
   m += NR_END;     int i;
   m -= nrl;     int lg=20;
       i=0;
       lg=strlen(t);
   /* allocate rows and set pointers to them */     for(i=0; i<= lg; i++) {
   m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));       (s[i] = t[i]);
   if (!m[nrl]) nrerror("allocation failure 2 in matrix()");       if (t[i]== '\\') s[i]='/';
   m[nrl] += NR_END;     }
   m[nrl] -= ncl;   }
     
   for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;   int nbocc(char *s, char occ)
     {
   /* return pointer to array of pointers to rows */     int i,j=0;
   return m;     int lg=20;
 }     i=0;
     lg=strlen(s);
 /****************** free_imatrix *************************/    for(i=0; i<= lg; i++) {
 void free_imatrix(m,nrl,nrh,ncl,nch)    if  (s[i] == occ ) j++;
       int **m;    }
       long nch,ncl,nrh,nrl;     return j;
      /* free an int matrix allocated by imatrix() */   }
 {   
   free((FREE_ARG) (m[nrl]+ncl-NR_END));   void cutv(char *u,char *v, char*t, char occ)
   free((FREE_ARG) (m+nrl-NR_END));   {
 }     int i,lg,j,p=0;
     i=0;
 /******************* matrix *******************************/    for(j=0; j<=strlen(t)-1; j++) {
 double **matrix(long nrl, long nrh, long ncl, long nch)      if((t[j]!= occ) && (t[j+1]== occ)) p=j+1;
 {    }
   long i, nrow=nrh-nrl+1, ncol=nch-ncl+1;  
   double **m;    lg=strlen(t);
     for(j=0; j<p; j++) {
   m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));      (u[j] = t[j]);
   if (!m) nrerror("allocation failure 1 in matrix()");    }
   m += NR_END;       u[p]='\0';
   m -= nrl;  
      for(j=0; j<= lg; j++) {
   m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));      if (j>=(p+1))(v[j-p-1] = t[j]);
   if (!m[nrl]) nrerror("allocation failure 2 in matrix()");    }
   m[nrl] += NR_END;  }
   m[nrl] -= ncl;  
   /********************** nrerror ********************/
   for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol;  
   return m;  void nrerror(char error_text[])
 }  {
     fprintf(stderr,"ERREUR ...\n");
 /*************************free matrix ************************/    fprintf(stderr,"%s\n",error_text);
 void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)    exit(1);
 {  }
   free((FREE_ARG)(m[nrl]+ncl-NR_END));  /*********************** vector *******************/
   free((FREE_ARG)(m+nrl-NR_END));  double *vector(int nl, int nh)
 }  {
     double *v;
 /******************* ma3x *******************************/    v=(double *) malloc((size_t)((nh-nl+1+NR_END)*sizeof(double)));
 double ***ma3x(long nrl, long nrh, long ncl, long nch, long nll, long nlh)    if (!v) nrerror("allocation failure in vector");
 {    return v-nl+NR_END;
   long i, j, nrow=nrh-nrl+1, ncol=nch-ncl+1, nlay=nlh-nll+1;  }
   double ***m;  
   /************************ free vector ******************/
   m=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double*)));  void free_vector(double*v, int nl, int nh)
   if (!m) nrerror("allocation failure 1 in matrix()");  {
   m += NR_END;    free((FREE_ARG)(v+nl-NR_END));
   m -= nrl;  }
   
   m[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));  /************************ivector *******************************/
   if (!m[nrl]) nrerror("allocation failure 2 in matrix()");  int *ivector(long nl,long nh)
   m[nrl] += NR_END;  {
   m[nrl] -= ncl;    int *v;
     v=(int *) malloc((size_t)((nh-nl+1+NR_END)*sizeof(int)));
   for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol;    if (!v) nrerror("allocation failure in ivector");
     return v-nl+NR_END;
   m[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*nlay+NR_END)*sizeof(double)));  }
   if (!m[nrl][ncl]) nrerror("allocation failure 3 in matrix()");  
   m[nrl][ncl] += NR_END;  /******************free ivector **************************/
   m[nrl][ncl] -= nll;  void free_ivector(int *v, long nl, long nh)
   for (j=ncl+1; j<=nch; j++)   {
     m[nrl][j]=m[nrl][j-1]+nlay;    free((FREE_ARG)(v+nl-NR_END));
     }
   for (i=nrl+1; i<=nrh; i++) {  
     m[i][ncl]=m[i-1l][ncl]+ncol*nlay;  /******************* imatrix *******************************/
     for (j=ncl+1; j<=nch; j++)   int **imatrix(long nrl, long nrh, long ncl, long nch)
       m[i][j]=m[i][j-1]+nlay;       /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
   }  {
   return m;    long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
 }    int **m;
    
 /*************************free ma3x ************************/    /* allocate pointers to rows */
 void free_ma3x(double ***m, long nrl, long nrh, long ncl, long nch,long nll, long nlh)    m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
 {    if (!m) nrerror("allocation failure 1 in matrix()");
   free((FREE_ARG)(m[nrl][ncl]+ nll-NR_END));    m += NR_END;
   free((FREE_ARG)(m[nrl]+ncl-NR_END));    m -= nrl;
   free((FREE_ARG)(m+nrl-NR_END));   
 }   
     /* allocate rows and set pointers to them */
 /***************** f1dim *************************/    m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
 extern int ncom;     if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
 extern double *pcom,*xicom;    m[nrl] += NR_END;
 extern double (*nrfunc)(double []);     m[nrl] -= ncl;
     
 double f1dim(double x)     for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
 {    
   int j;     /* return pointer to array of pointers to rows */
   double f;    return m;
   double *xt;   }
    
   xt=vector(1,ncom);   /****************** free_imatrix *************************/
   for (j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];   void free_imatrix(m,nrl,nrh,ncl,nch)
   f=(*nrfunc)(xt);         int **m;
   free_vector(xt,1,ncom);         long nch,ncl,nrh,nrl;
   return f;        /* free an int matrix allocated by imatrix() */
 }   {
     free((FREE_ARG) (m[nrl]+ncl-NR_END));
 /*****************brent *************************/    free((FREE_ARG) (m+nrl-NR_END));
 double brent(double ax, double bx, double cx, double (*f)(double), double tol,  double *xmin)   }
 {   
   int iter;   /******************* matrix *******************************/
   double a,b,d,etemp;  double **matrix(long nrl, long nrh, long ncl, long nch)
   double fu,fv,fw,fx;  {
   double ftemp;    long i, nrow=nrh-nrl+1, ncol=nch-ncl+1;
   double p,q,r,tol1,tol2,u,v,w,x,xm;     double **m;
   double e=0.0;   
      m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
   a=(ax < cx ? ax : cx);     if (!m) nrerror("allocation failure 1 in matrix()");
   b=(ax > cx ? ax : cx);     m += NR_END;
   x=w=v=bx;     m -= nrl;
   fw=fv=fx=(*f)(x);   
   for (iter=1;iter<=ITMAX;iter++) {     m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
     xm=0.5*(a+b);     if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
     tol2=2.0*(tol1=tol*fabs(x)+ZEPS);     m[nrl] += NR_END;
     /*          if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret)))*/    m[nrl] -= ncl;
     printf(".");fflush(stdout);  
 #ifdef DEBUG    for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol;
     printf("br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);    return m;
     /*          if ((fabs(x-xm) <= (tol2-0.5*(b-a)))||(2.0*fabs(fu-ftemp) <= ftol*1.e-2*(fabs(fu)+fabs(ftemp)))) { */  }
 #endif  
     if (fabs(x-xm) <= (tol2-0.5*(b-a))){   /*************************free matrix ************************/
       *xmin=x;   void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
       return fx;   {
     }     free((FREE_ARG)(m[nrl]+ncl-NR_END));
     ftemp=fu;    free((FREE_ARG)(m+nrl-NR_END));
     if (fabs(e) > tol1) {   }
       r=(x-w)*(fx-fv);   
       q=(x-v)*(fx-fw);   /******************* ma3x *******************************/
       p=(x-v)*q-(x-w)*r;   double ***ma3x(long nrl, long nrh, long ncl, long nch, long nll, long nlh)
       q=2.0*(q-r);   {
       if (q > 0.0) p = -p;     long i, j, nrow=nrh-nrl+1, ncol=nch-ncl+1, nlay=nlh-nll+1;
       q=fabs(q);     double ***m;
       etemp=e;   
       e=d;     m=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
       if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))     if (!m) nrerror("allocation failure 1 in matrix()");
         d=CGOLD*(e=(x >= xm ? a-x : b-x));     m += NR_END;
       else {     m -= nrl;
         d=p/q;   
         u=x+d;     m[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
         if (u-a < tol2 || b-u < tol2)     if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
           d=SIGN(tol1,xm-x);     m[nrl] += NR_END;
       }     m[nrl] -= ncl;
     } else {   
       d=CGOLD*(e=(x >= xm ? a-x : b-x));     for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol;
     }   
     u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));     m[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*nlay+NR_END)*sizeof(double)));
     fu=(*f)(u);     if (!m[nrl][ncl]) nrerror("allocation failure 3 in matrix()");
     if (fu <= fx) {     m[nrl][ncl] += NR_END;
       if (u >= x) a=x; else b=x;     m[nrl][ncl] -= nll;
       SHFT(v,w,x,u)     for (j=ncl+1; j<=nch; j++)
         SHFT(fv,fw,fx,fu)       m[nrl][j]=m[nrl][j-1]+nlay;
         } else {    
           if (u < x) a=u; else b=u;     for (i=nrl+1; i<=nrh; i++) {
           if (fu <= fw || w == x) {       m[i][ncl]=m[i-1l][ncl]+ncol*nlay;
             v=w;       for (j=ncl+1; j<=nch; j++)
             w=u;         m[i][j]=m[i][j-1]+nlay;
             fv=fw;     }
             fw=fu;     return m;
           } else if (fu <= fv || v == x || v == w) {   }
             v=u;   
             fv=fu;   /*************************free ma3x ************************/
           }   void free_ma3x(double ***m, long nrl, long nrh, long ncl, long nch,long nll, long nlh)
         }   {
   }     free((FREE_ARG)(m[nrl][ncl]+ nll-NR_END));
   nrerror("Too many iterations in brent");     free((FREE_ARG)(m[nrl]+ncl-NR_END));
   *xmin=x;     free((FREE_ARG)(m+nrl-NR_END));
   return fx;   }
 }   
   /***************** f1dim *************************/
 /****************** mnbrak ***********************/  extern int ncom;
   extern double *pcom,*xicom;
 void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc,   extern double (*nrfunc)(double []);
             double (*func)(double))    
 {   double f1dim(double x)
   double ulim,u,r,q, dum;  {
   double fu;     int j;
      double f;
   *fa=(*func)(*ax);     double *xt;
   *fb=(*func)(*bx);    
   if (*fb > *fa) {     xt=vector(1,ncom);
     SHFT(dum,*ax,*bx,dum)     for (j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];
       SHFT(dum,*fb,*fa,dum)     f=(*nrfunc)(xt);
       }     free_vector(xt,1,ncom);
   *cx=(*bx)+GOLD*(*bx-*ax);     return f;
   *fc=(*func)(*cx);   }
   while (*fb > *fc) {   
     r=(*bx-*ax)*(*fb-*fc);   /*****************brent *************************/
     q=(*bx-*cx)*(*fb-*fa);   double brent(double ax, double bx, double cx, double (*f)(double), double tol,  double *xmin)
     u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/   {
       (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));     int iter;
     ulim=(*bx)+GLIMIT*(*cx-*bx);     double a,b,d,etemp;
     if ((*bx-u)*(u-*cx) > 0.0) {     double fu,fv,fw,fx;
       fu=(*func)(u);     double ftemp;
     } else if ((*cx-u)*(u-ulim) > 0.0) {     double p,q,r,tol1,tol2,u,v,w,x,xm;
       fu=(*func)(u);     double e=0.0;
       if (fu < *fc) {    
         SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))     a=(ax < cx ? ax : cx);
           SHFT(*fb,*fc,fu,(*func)(u))     b=(ax > cx ? ax : cx);
           }     x=w=v=bx;
     } else if ((u-ulim)*(ulim-*cx) >= 0.0) {     fw=fv=fx=(*f)(x);
       u=ulim;     for (iter=1;iter<=ITMAX;iter++) {
       fu=(*func)(u);       xm=0.5*(a+b);
     } else {       tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
       u=(*cx)+GOLD*(*cx-*bx);       /*          if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret)))*/
       fu=(*func)(u);       printf(".");fflush(stdout);
     }   #ifdef DEBUG
     SHFT(*ax,*bx,*cx,u)       printf("br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);
       SHFT(*fa,*fb,*fc,fu)       /*          if ((fabs(x-xm) <= (tol2-0.5*(b-a)))||(2.0*fabs(fu-ftemp) <= ftol*1.e-2*(fabs(fu)+fabs(ftemp)))) { */
       }   #endif
 }       if (fabs(x-xm) <= (tol2-0.5*(b-a))){
         *xmin=x;
 /*************** linmin ************************/        return fx;
       }
 int ncom;       ftemp=fu;
 double *pcom,*xicom;      if (fabs(e) > tol1) {
 double (*nrfunc)(double []);         r=(x-w)*(fx-fv);
          q=(x-v)*(fx-fw);
 void linmin(double p[], double xi[], int n, double *fret,double (*func)(double []))         p=(x-v)*q-(x-w)*r;
 {         q=2.0*(q-r);
   double brent(double ax, double bx, double cx,         if (q > 0.0) p = -p;
                double (*f)(double), double tol, double *xmin);         q=fabs(q);
   double f1dim(double x);         etemp=e;
   void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb,         e=d;
               double *fc, double (*func)(double));         if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
   int j;           d=CGOLD*(e=(x >= xm ? a-x : b-x));
   double xx,xmin,bx,ax;         else {
   double fx,fb,fa;          d=p/q;
            u=x+d;
   ncom=n;           if (u-a < tol2 || b-u < tol2)
   pcom=vector(1,n);             d=SIGN(tol1,xm-x);
   xicom=vector(1,n);         }
   nrfunc=func;       } else {
   for (j=1;j<=n;j++) {         d=CGOLD*(e=(x >= xm ? a-x : b-x));
     pcom[j]=p[j];       }
     xicom[j]=xi[j];       u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
   }       fu=(*f)(u);
   ax=0.0;       if (fu <= fx) {
   xx=1.0;         if (u >= x) a=x; else b=x;
   mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);         SHFT(v,w,x,u)
   *fret=brent(ax,xx,bx,f1dim,TOL,&xmin);           SHFT(fv,fw,fx,fu)
 #ifdef DEBUG          } else {
   printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);            if (u < x) a=u; else b=u;
 #endif            if (fu <= fw || w == x) {
   for (j=1;j<=n;j++) {               v=w;
     xi[j] *= xmin;               w=u;
     p[j] += xi[j];               fv=fw;
   }               fw=fu;
   free_vector(xicom,1,n);             } else if (fu <= fv || v == x || v == w) {
   free_vector(pcom,1,n);               v=u;
 }               fv=fu;
             }
 /*************** powell ************************/          }
 void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret,     }
             double (*func)(double []))     nrerror("Too many iterations in brent");
     *xmin=x;
 {     return fx;
   }
   
   void linmin(double p[], double xi[], int n, double *fret,   /****************** mnbrak ***********************/
               double (*func)(double []));   
   int i,ibig,j;   void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc,
   double del,t,*pt,*ptt,*xit;              double (*func)(double))
   double fp,fptt;  {
   double *xits;    double ulim,u,r,q, dum;
   pt=vector(1,n);     double fu;
   ptt=vector(1,n);    
   xit=vector(1,n);     *fa=(*func)(*ax);
   xits=vector(1,n);     *fb=(*func)(*bx);
   *fret=(*func)(p);     if (*fb > *fa) {
   for (j=1;j<=n;j++) pt[j]=p[j];       SHFT(dum,*ax,*bx,dum)
   for (*iter=1;;++(*iter)) {         SHFT(dum,*fb,*fa,dum)
     fp=(*fret);         }
     ibig=0;     *cx=(*bx)+GOLD*(*bx-*ax);
     del=0.0;     *fc=(*func)(*cx);
     printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret);    while (*fb > *fc) {
     for (i=1;i<=n;i++)       r=(*bx-*ax)*(*fb-*fc);
       printf(" %d %.12f",i, p[i]);      q=(*bx-*cx)*(*fb-*fa);
     printf("\n");      u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
     for (i=1;i<=n;i++) {         (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));
       for (j=1;j<=n;j++) xit[j]=xi[j][i];       ulim=(*bx)+GLIMIT*(*cx-*bx);
       fptt=(*fret);       if ((*bx-u)*(u-*cx) > 0.0) {
 #ifdef DEBUG        fu=(*func)(u);
       printf("fret=%lf \n",*fret);      } else if ((*cx-u)*(u-ulim) > 0.0) {
 #endif        fu=(*func)(u);
       printf("%d",i);fflush(stdout);        if (fu < *fc) {
       linmin(p,xit,n,fret,func);           SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
       if (fabs(fptt-(*fret)) > del) {             SHFT(*fb,*fc,fu,(*func)(u))
         del=fabs(fptt-(*fret));             }
         ibig=i;       } else if ((u-ulim)*(ulim-*cx) >= 0.0) {
       }         u=ulim;
 #ifdef DEBUG        fu=(*func)(u);
       printf("%d %.12e",i,(*fret));      } else {
       for (j=1;j<=n;j++) {        u=(*cx)+GOLD*(*cx-*bx);
         xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);        fu=(*func)(u);
         printf(" x(%d)=%.12e",j,xit[j]);      }
       }      SHFT(*ax,*bx,*cx,u)
       for(j=1;j<=n;j++)         SHFT(*fa,*fb,*fc,fu)
         printf(" p=%.12e",p[j]);        }
       printf("\n");  }
 #endif  
     }   /*************** linmin ************************/
     if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {  
 #ifdef DEBUG  int ncom;
       int k[2],l;  double *pcom,*xicom;
       k[0]=1;  double (*nrfunc)(double []);
       k[1]=-1;   
       printf("Max: %.12e",(*func)(p));  void linmin(double p[], double xi[], int n, double *fret,double (*func)(double []))
       for (j=1;j<=n;j++)   {
         printf(" %.12e",p[j]);    double brent(double ax, double bx, double cx,
       printf("\n");                 double (*f)(double), double tol, double *xmin);
       for(l=0;l<=1;l++) {    double f1dim(double x);
         for (j=1;j<=n;j++) {    void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb,
           ptt[j]=p[j]+(p[j]-pt[j])*k[l];                double *fc, double (*func)(double));
           printf("l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);    int j;
         }    double xx,xmin,bx,ax;
         printf("func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));    double fx,fb,fa;
       }   
 #endif    ncom=n;
     pcom=vector(1,n);
     xicom=vector(1,n);
       free_vector(xit,1,n);     nrfunc=func;
       free_vector(xits,1,n);     for (j=1;j<=n;j++) {
       free_vector(ptt,1,n);       pcom[j]=p[j];
       free_vector(pt,1,n);       xicom[j]=xi[j];
       return;     }
     }     ax=0.0;
     if (*iter == ITMAX) nrerror("powell exceeding maximum iterations.");     xx=1.0;
     for (j=1;j<=n;j++) {     mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);
       ptt[j]=2.0*p[j]-pt[j];     *fret=brent(ax,xx,bx,f1dim,TOL,&xmin);
       xit[j]=p[j]-pt[j];   #ifdef DEBUG
       pt[j]=p[j];     printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
     }   #endif
     fptt=(*func)(ptt);     for (j=1;j<=n;j++) {
     if (fptt < fp) {       xi[j] *= xmin;
       t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt);       p[j] += xi[j];
       if (t < 0.0) {     }
         linmin(p,xit,n,fret,func);     free_vector(xicom,1,n);
         for (j=1;j<=n;j++) {     free_vector(pcom,1,n);
           xi[j][ibig]=xi[j][n];   }
           xi[j][n]=xit[j];   
         }  /*************** powell ************************/
 #ifdef DEBUG  void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret,
         printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);              double (*func)(double []))
         for(j=1;j<=n;j++)  {
           printf(" %.12e",xit[j]);    void linmin(double p[], double xi[], int n, double *fret,
         printf("\n");                double (*func)(double []));
 #endif    int i,ibig,j;
       }     double del,t,*pt,*ptt,*xit;
     }     double fp,fptt;
   }     double *xits;
 }     pt=vector(1,n);
     ptt=vector(1,n);
 /**** Prevalence limit ****************/    xit=vector(1,n);
     xits=vector(1,n);
 double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl)    *fret=(*func)(p);
 {    for (j=1;j<=n;j++) pt[j]=p[j];
   /* Computes the prevalence limit in each live state at age x by left multiplying the unit    for (*iter=1;;++(*iter)) {
      matrix by transitions matrix until convergence is reached */      fp=(*fret);
       ibig=0;
   int i, ii,j,k;      del=0.0;
   double min, max, maxmin, maxmax,sumnew=0.;      printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret);
   double **matprod2();      for (i=1;i<=n;i++)
   double **out, cov[NCOVMAX], **pmij();        printf(" %d %.12f",i, p[i]);
   double **newm;      printf("\n");
   double agefin, delaymax=50 ; /* Max number of years to converge */      for (i=1;i<=n;i++) {
         for (j=1;j<=n;j++) xit[j]=xi[j][i];
   for (ii=1;ii<=nlstate+ndeath;ii++)        fptt=(*fret);
     for (j=1;j<=nlstate+ndeath;j++){  #ifdef DEBUG
       oldm[ii][j]=(ii==j ? 1.0 : 0.0);        printf("fret=%lf \n",*fret);
     }  #endif
   /* Even if hstepm = 1, at least one multiplication by the unit matrix */        printf("%d",i);fflush(stdout);
   for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){        linmin(p,xit,n,fret,func);
     newm=savm;        if (fabs(fptt-(*fret)) > del) {
     /* Covariates have to be included here again */          del=fabs(fptt-(*fret));
     cov[1]=1.;          ibig=i;
     cov[2]=agefin;        }
     out=matprod2(newm, pmij(pmmij,cov,ncov,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);  #ifdef DEBUG
 /*    printf("age=%f agefin=%f po=%f pn=%f\n",age,agefin,oldm[1][1],newm[1][1]);*/        printf("%d %.12e",i,(*fret));
             for (j=1;j<=n;j++) {
     savm=oldm;          xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
     oldm=newm;          printf(" x(%d)=%.12e",j,xit[j]);
     maxmax=0.;        }
     for(j=1;j<=nlstate;j++){        for(j=1;j<=n;j++)
       min=1.;          printf(" p=%.12e",p[j]);
       max=0.;        printf("\n");
       for(i=1; i<=nlstate; i++) {  #endif
         sumnew=0;      }
         for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];      if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
         prlim[i][j]= newm[i][j]/(1-sumnew);  #ifdef DEBUG
         max=FMAX(max,prlim[i][j]);        int k[2],l;
         min=FMIN(min,prlim[i][j]);        k[0]=1;
       }        k[1]=-1;
       maxmin=max-min;        printf("Max: %.12e",(*func)(p));
       maxmax=FMAX(maxmax,maxmin);        for (j=1;j<=n;j++)
     }          printf(" %.12e",p[j]);
     if(maxmax < ftolpl){        printf("\n");
       return prlim;        for(l=0;l<=1;l++) {
     }          for (j=1;j<=n;j++) {
   }            ptt[j]=p[j]+(p[j]-pt[j])*k[l];
 }            printf("l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
           }
 /*************** transition probabilities **********/           printf("func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));
         }
 double **pmij(double **ps, double *cov, int ncov, double *x, int nlstate )  #endif
 {  
   double s1, s2;  
   /*double t34;*/        free_vector(xit,1,n);
   int i,j,j1, nc, ii, jj;        free_vector(xits,1,n);
         free_vector(ptt,1,n);
     for(i=1; i<= nlstate; i++){        free_vector(pt,1,n);
     for(j=1; j<i;j++){        return;
       for (nc=1, s2=0.;nc <=ncov; nc++){      }
         /*s2 += param[i][j][nc]*cov[nc];*/      if (*iter == ITMAX) nrerror("powell exceeding maximum iterations.");
         s2 += x[(i-1)*nlstate*ncov+(j-1)*ncov+nc+(i-1)*(ndeath-1)*ncov]*cov[nc];      for (j=1;j<=n;j++) {
         /*printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2);*/        ptt[j]=2.0*p[j]-pt[j];
       }        xit[j]=p[j]-pt[j];
       ps[i][j]=s2;        pt[j]=p[j];
       /*printf("s1=%.17e, s2=%.17e\n",s1,s2);*/      }
     }      fptt=(*func)(ptt);
     for(j=i+1; j<=nlstate+ndeath;j++){      if (fptt < fp) {
       for (nc=1, s2=0.;nc <=ncov; nc++){        t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt);
         s2 += x[(i-1)*nlstate*ncov+(j-2)*ncov+nc+(i-1)*(ndeath-1)*ncov]*cov[nc];        if (t < 0.0) {
         /*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/          linmin(p,xit,n,fret,func);
       }          for (j=1;j<=n;j++) {
       ps[i][j]=s2;            xi[j][ibig]=xi[j][n];
     }            xi[j][n]=xit[j];
   }          }
   for(i=1; i<= nlstate; i++){  #ifdef DEBUG
      s1=0;          printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
     for(j=1; j<i; j++)          for(j=1;j<=n;j++)
       s1+=exp(ps[i][j]);            printf(" %.12e",xit[j]);
     for(j=i+1; j<=nlstate+ndeath; j++)          printf("\n");
       s1+=exp(ps[i][j]);  #endif
     ps[i][i]=1./(s1+1.);        }
     for(j=1; j<i; j++)      }
       ps[i][j]= exp(ps[i][j])*ps[i][i];    }
     for(j=i+1; j<=nlstate+ndeath; j++)  }
       ps[i][j]= exp(ps[i][j])*ps[i][i];  
     /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */  /**** Prevalence limit ****************/
   } /* end i */  
   double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij)
   for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){  {
     for(jj=1; jj<= nlstate+ndeath; jj++){    /* Computes the prevalence limit in each live state at age x by left multiplying the unit
       ps[ii][jj]=0;       matrix by transitions matrix until convergence is reached */
       ps[ii][ii]=1;  
     }    int i, ii,j,k;
   }    double min, max, maxmin, maxmax,sumnew=0.;
     double **matprod2();
   /*   for(ii=1; ii<= nlstate+ndeath; ii++){    double **out, cov[NCOVMAX], **pmij();
     for(jj=1; jj<= nlstate+ndeath; jj++){    double **newm;
      printf("%lf ",ps[ii][jj]);    double agefin, delaymax=50 ; /* Max number of years to converge */
    }  
     printf("\n ");    for (ii=1;ii<=nlstate+ndeath;ii++)
     }      for (j=1;j<=nlstate+ndeath;j++){
     printf("\n ");printf("%lf ",cov[2]);*/        oldm[ii][j]=(ii==j ? 1.0 : 0.0);
 /*      }
   for(i=1; i<= npar; i++) printf("%f ",x[i]);  
   goto end;*/     cov[1]=1.;
     return ps;   
 }   /* Even if hstepm = 1, at least one multiplication by the unit matrix */
     for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
 /**************** Product of 2 matrices ******************/      newm=savm;
       /* Covariates have to be included here again */
 double **matprod2(double **out, double **in,long nrl, long nrh, long ncl, long nch, long ncolol, long ncoloh, double **b)       cov[2]=agefin;
 {   
   /* Computes the matric product of in(1,nrh-nrl+1)(1,nch-ncl+1) times        for (k=1; k<=cptcovn;k++) {
      b(1,nch-ncl+1)(1,ncoloh-ncolol+1) into out(...) */          cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
   /* in, b, out are matrice of pointers which should have been initialized           /*      printf("ij=%d k=%d Tvar[k]=%d nbcode=%d cov=%lf codtab[ij][Tvar[k]]=%d \n",ij,k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], codtab[ij][Tvar[k]]);*/
      before: only the contents of out is modified. The function returns        }
      a pointer to pointers identical to out */        for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
   long i, j, k;        for (k=1; k<=cptcovprod;k++)
   for(i=nrl; i<= nrh; i++)          cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
     for(k=ncolol; k<=ncoloh; k++)  
       for(j=ncl,out[i][k]=0.; j<=nch; j++)        /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
         out[i][k] +=in[i][j]*b[j][k];        /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
         /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
   return out;      out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
 }  
       savm=oldm;
       oldm=newm;
 /************* Higher Matrix Product ***************/      maxmax=0.;
       for(j=1;j<=nlstate;j++){
 double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm )        min=1.;
 {        max=0.;
   /* Computes the transition matrix starting at age 'age' over 'nhstepm*hstepm*stepm' month         for(i=1; i<=nlstate; i++) {
      duration (i.e. until          sumnew=0;
      age (in years)  age+nhstepm*stepm/12) by multiplying nhstepm*hstepm matrices.           for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];
      Output is stored in matrix po[i][j][h] for h every 'hstepm' step           prlim[i][j]= newm[i][j]/(1-sumnew);
      (typically every 2 years instead of every month which is too big).          max=FMAX(max,prlim[i][j]);
      Model is determined by parameters x and covariates have to be           min=FMIN(min,prlim[i][j]);
      included manually here.         }
         maxmin=max-min;
      */        maxmax=FMAX(maxmax,maxmin);
       }
   int i, j, d, h;      if(maxmax < ftolpl){
   double **out, cov[NCOVMAX];        return prlim;
   double **newm;      }
     }
   /* Hstepm could be zero and should return the unit matrix */  }
   for (i=1;i<=nlstate+ndeath;i++)  
     for (j=1;j<=nlstate+ndeath;j++){  /*************** transition probabilities ***************/
       oldm[i][j]=(i==j ? 1.0 : 0.0);  
       po[i][j][0]=(i==j ? 1.0 : 0.0);  double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
     }  {
   /* Even if hstepm = 1, at least one multiplication by the unit matrix */    double s1, s2;
   for(h=1; h <=nhstepm; h++){    /*double t34;*/
     for(d=1; d <=hstepm; d++){    int i,j,j1, nc, ii, jj;
       newm=savm;  
       /* Covariates have to be included here again */      for(i=1; i<= nlstate; i++){
       cov[1]=1.;      for(j=1; j<i;j++){
       cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;        for (nc=1, s2=0.;nc <=ncovmodel; nc++){
       /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/          /*s2 += param[i][j][nc]*cov[nc];*/
       out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,           s2 += x[(i-1)*nlstate*ncovmodel+(j-1)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
                    pmij(pmmij,cov,ncov,x,nlstate));          /*printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2);*/
       savm=oldm;        }
       oldm=newm;        ps[i][j]=s2;
     }        /*printf("s1=%.17e, s2=%.17e\n",s1,s2);*/
     for(i=1; i<=nlstate+ndeath; i++)      }
       for(j=1;j<=nlstate+ndeath;j++) {      for(j=i+1; j<=nlstate+ndeath;j++){
         po[i][j][h]=newm[i][j];        for (nc=1, s2=0.;nc <=ncovmodel; nc++){
         /*printf("i=%d j=%d h=%d po[i][j][h]=%f ",i,j,h,po[i][j][h]);          s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
          */          /*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
       }        }
   } /* end h */        ps[i][j]=s2;
   return po;      }
 }    }
       /*ps[3][2]=1;*/
   
 /*************** log-likelihood *************/    for(i=1; i<= nlstate; i++){
 double func( double *x)       s1=0;
 {      for(j=1; j<i; j++)
   int i, ii, j, k, mi, d;        s1+=exp(ps[i][j]);
   double l, ll[NLSTATEMAX], cov[NCOVMAX];      for(j=i+1; j<=nlstate+ndeath; j++)
   double **out;        s1+=exp(ps[i][j]);
   double sw; /* Sum of weights */      ps[i][i]=1./(s1+1.);
   double lli; /* Individual log likelihood */      for(j=1; j<i; j++)
   long ipmx;        ps[i][j]= exp(ps[i][j])*ps[i][i];
   /*extern weight */      for(j=i+1; j<=nlstate+ndeath; j++)
   /* We are differentiating ll according to initial status */        ps[i][j]= exp(ps[i][j])*ps[i][i];
   /*  for (i=1;i<=npar;i++) printf("%f ", x[i]);*/      /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
   /*for(i=1;i<imx;i++)     } /* end i */
 printf(" %d\n",s[4][i]);  
   */    for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
       for(jj=1; jj<= nlstate+ndeath; jj++){
   for(k=1; k<=nlstate; k++) ll[k]=0.;        ps[ii][jj]=0;
   for (i=1,ipmx=0, sw=0.; i<=imx; i++){        ps[ii][ii]=1;
        for(mi=1; mi<= wav[i]-1; mi++){      }
       for (ii=1;ii<=nlstate+ndeath;ii++)    }
         for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0);  
             for(d=0; d<dh[mi][i]; d++){  
         newm=savm;    /*   for(ii=1; ii<= nlstate+ndeath; ii++){
           cov[1]=1.;      for(jj=1; jj<= nlstate+ndeath; jj++){
           cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;       printf("%lf ",ps[ii][jj]);
             out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,     }
                        1,nlstate+ndeath,pmij(pmmij,cov,ncov,x,nlstate));      printf("\n ");
           savm=oldm;      }
           oldm=newm;      printf("\n ");printf("%lf ",cov[2]);*/
   /*
     for(i=1; i<= npar; i++) printf("%f ",x[i]);
       } /* end mult */    goto end;*/
          return ps;
       lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);  }
       /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/  
       ipmx +=1;  /**************** Product of 2 matrices ******************/
       sw += weight[i];  
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;  double **matprod2(double **out, double **in,long nrl, long nrh, long ncl, long nch, long ncolol, long ncoloh, double **b)
     } /* end of wave */  {
   } /* end of individual */    /* 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(...) */
   for(k=1,l=0.; k<=nlstate; k++) l += ll[k];    /* in, b, out are matrice of pointers which should have been initialized
   /* printf("l1=%f l2=%f ",ll[1],ll[2]); */       before: only the contents of out is modified. The function returns
   l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */       a pointer to pointers identical to out */
   return -l;    long i, j, k;
 }    for(i=nrl; i<= nrh; i++)
       for(k=ncolol; k<=ncoloh; k++)
         for(j=ncl,out[i][k]=0.; j<=nch; j++)
 /*********** Maximum Likelihood Estimation ***************/          out[i][k] +=in[i][j]*b[j][k];
   
 void mlikeli(FILE *ficres,double p[], int npar, int ncov, int nlstate, double ftol, double (*func)(double []))    return out;
 {  }
   int i,j, iter;  
   double **xi,*delti;  
   double fret;  /************* Higher Matrix Product ***************/
   xi=matrix(1,npar,1,npar);  
   for (i=1;i<=npar;i++)  double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij )
     for (j=1;j<=npar;j++)  {
       xi[i][j]=(i==j ? 1.0 : 0.0);    /* Computes the transition matrix starting at age 'age' over 'nhstepm*hstepm*stepm' month
   printf("Powell\n");       duration (i.e. until
   powell(p,xi,npar,ftol,&iter,&fret,func);       age (in years)  age+nhstepm*stepm/12) by multiplying nhstepm*hstepm matrices.
        Output is stored in matrix po[i][j][h] for h every 'hstepm' step
    printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));       (typically every 2 years instead of every month which is too big).
   fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f ",iter,func(p));       Model is determined by parameters x and covariates have to be
        included manually here.
 }  
        */
 /**** Computes Hessian and covariance matrix ***/  
 void hesscov(double **matcov, double p[], int npar, double delti[], double ftolhess, double (*func)(double []))    int i, j, d, h, k;
 {    double **out, cov[NCOVMAX];
   double  **a,**y,*x,pd;    double **newm;
   double **hess;  
   int i, j,jk;    /* Hstepm could be zero and should return the unit matrix */
   int *indx;    for (i=1;i<=nlstate+ndeath;i++)
       for (j=1;j<=nlstate+ndeath;j++){
   double hessii(double p[], double delta, int theta, double delti[]);        oldm[i][j]=(i==j ? 1.0 : 0.0);
   double hessij(double p[], double delti[], int i, int j);        po[i][j][0]=(i==j ? 1.0 : 0.0);
   void lubksb(double **a, int npar, int *indx, double b[]) ;      }
   void ludcmp(double **a, int npar, int *indx, double *d) ;    /* Even if hstepm = 1, at least one multiplication by the unit matrix */
     for(h=1; h <=nhstepm; h++){
       for(d=1; d <=hstepm; d++){
   hess=matrix(1,npar,1,npar);        newm=savm;
         /* Covariates have to be included here again */
   printf("\nCalculation of the hessian matrix. Wait...\n");        cov[1]=1.;
   for (i=1;i<=npar;i++){        cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
     printf("%d",i);fflush(stdout);        for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
     hess[i][i]=hessii(p,ftolhess,i,delti);        for (k=1; k<=cptcovage;k++)
     /*printf(" %f ",p[i]);*/          cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
   }        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]]];
   for (i=1;i<=npar;i++) {  
     for (j=1;j<=npar;j++)  {  
       if (j>i) {         /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
         printf(".%d%d",i,j);fflush(stdout);        /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
         hess[i][j]=hessij(p,delti,i,j);        out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,
         hess[j][i]=hess[i][j];                     pmij(pmmij,cov,ncovmodel,x,nlstate));
       }        savm=oldm;
     }        oldm=newm;
   }      }
   printf("\n");      for(i=1; i<=nlstate+ndeath; i++)
         for(j=1;j<=nlstate+ndeath;j++) {
   printf("\nInverting the hessian to get the covariance matrix. Wait...\n");          po[i][j][h]=newm[i][j];
             /*printf("i=%d j=%d h=%d po[i][j][h]=%f ",i,j,h,po[i][j][h]);
   a=matrix(1,npar,1,npar);           */
   y=matrix(1,npar,1,npar);        }
   x=vector(1,npar);    } /* end h */
   indx=ivector(1,npar);    return po;
   for (i=1;i<=npar;i++)  }
     for (j=1;j<=npar;j++) a[i][j]=hess[i][j];  
   ludcmp(a,npar,indx,&pd);  
   /*************** log-likelihood *************/
   for (j=1;j<=npar;j++) {  double func( double *x)
     for (i=1;i<=npar;i++) x[i]=0;  {
     x[j]=1;    int i, ii, j, k, mi, d, kk;
     lubksb(a,npar,indx,x);    double l, ll[NLSTATEMAX], cov[NCOVMAX];
     for (i=1;i<=npar;i++){     double **out;
       matcov[i][j]=x[i];    double sw; /* Sum of weights */
     }    double lli; /* Individual log likelihood */
   }    long ipmx;
     /*extern weight */
   printf("\n#Hessian matrix#\n");    /* We are differentiating ll according to initial status */
   for (i=1;i<=npar;i++) {     /*  for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
     for (j=1;j<=npar;j++) {     /*for(i=1;i<imx;i++)
       printf("%.3e ",hess[i][j]);      printf(" %d\n",s[4][i]);
     }    */
     printf("\n");    cov[1]=1.;
   }  
     for(k=1; k<=nlstate; k++) ll[k]=0.;
   /* Recompute Inverse */    for (i=1,ipmx=0, sw=0.; i<=imx; i++){
   for (i=1;i<=npar;i++)      for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
     for (j=1;j<=npar;j++) a[i][j]=matcov[i][j];      for(mi=1; mi<= wav[i]-1; mi++){
   ludcmp(a,npar,indx,&pd);        for (ii=1;ii<=nlstate+ndeath;ii++)
           for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0);
   /*  printf("\n#Hessian matrix recomputed#\n");        for(d=0; d<dh[mi][i]; d++){
           newm=savm;
   for (j=1;j<=npar;j++) {          cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
     for (i=1;i<=npar;i++) x[i]=0;          for (kk=1; kk<=cptcovage;kk++) {
     x[j]=1;            cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
     lubksb(a,npar,indx,x);          }
     for (i=1;i<=npar;i++){          
       y[i][j]=x[i];          out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
       printf("%.3e ",y[i][j]);                       1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
     }          savm=oldm;
     printf("\n");          oldm=newm;
   }         
   */         
         } /* end mult */
   free_matrix(a,1,npar,1,npar);       
   free_matrix(y,1,npar,1,npar);        lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);
   free_vector(x,1,npar);        /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/
   free_ivector(indx,1,npar);        ipmx +=1;
   free_matrix(hess,1,npar,1,npar);        sw += weight[i];
         ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
       } /* end of wave */
 }    } /* end of individual */
   
 /*************** hessian matrix ****************/    for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
 double hessii( double x[], double delta, int theta, double delti[])    /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
 {    l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
   int i;    return -l;
   int l=1, lmax=20;  }
   double k1,k2;  
   double p2[NPARMAX+1];  
   double res;  /*********** Maximum Likelihood Estimation ***************/
   double delt, delts, nkhi=10.,nkhif=1., khi=1.e-4;  
   double fx;  void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double []))
   int k=0,kmax=10;  {
   double l1;    int i,j, iter;
     double **xi,*delti;
   fx=func(x);    double fret;
   for (i=1;i<=npar;i++) p2[i]=x[i];    xi=matrix(1,npar,1,npar);
   for(l=0 ; l <=lmax; l++){    for (i=1;i<=npar;i++)
     l1=pow(10,l);      for (j=1;j<=npar;j++)
     delts=delt;        xi[i][j]=(i==j ? 1.0 : 0.0);
     for(k=1 ; k <kmax; k=k+1){    printf("Powell\n");
       delt = delta*(l1*k);    powell(p,xi,npar,ftol,&iter,&fret,func);
       p2[theta]=x[theta] +delt;  
       k1=func(p2)-fx;     printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
       p2[theta]=x[theta]-delt;    fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
       k2=func(p2)-fx;  
       /*res= (k1-2.0*fx+k2)/delt/delt; */  }
       res= (k1+k2)/delt/delt/2.; /* Divided by because L and not 2*L */  
         /**** Computes Hessian and covariance matrix ***/
 #ifdef DEBUG  void hesscov(double **matcov, double p[], int npar, double delti[], double ftolhess, double (*func)(double []))
       printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx);  {
 #endif    double  **a,**y,*x,pd;
       /*if(fabs(k1-2.0*fx+k2) <1.e-13){ */    double **hess;
       if((k1 <khi/nkhi/2.) || (k2 <khi/nkhi/2.)){    int i, j,jk;
         k=kmax;    int *indx;
       }  
       else if((k1 >khi/nkhif) || (k2 >khi/nkhif)){ /* Keeps lastvalue before 3.84/2 KHI2 5% 1d.f. */    double hessii(double p[], double delta, int theta, double delti[]);
         k=kmax; l=lmax*10.;    double hessij(double p[], double delti[], int i, int j);
       }    void lubksb(double **a, int npar, int *indx, double b[]) ;
       else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){     void ludcmp(double **a, int npar, int *indx, double *d) ;
         delts=delt;  
       }    hess=matrix(1,npar,1,npar);
     }  
   }    printf("\nCalculation of the hessian matrix. Wait...\n");
   delti[theta]=delts;    for (i=1;i<=npar;i++){
   return res;      printf("%d",i);fflush(stdout);
         hess[i][i]=hessii(p,ftolhess,i,delti);
 }      /*printf(" %f ",p[i]);*/
       /*printf(" %lf ",hess[i][i]);*/
 double hessij( double x[], double delti[], int thetai,int thetaj)    }
 {   
   int i;    for (i=1;i<=npar;i++) {
   int l=1, l1, lmax=20;      for (j=1;j<=npar;j++)  {
   double k1,k2,k3,k4,res,fx;        if (j>i) {
   double p2[NPARMAX+1];          printf(".%d%d",i,j);fflush(stdout);
   int k;          hess[i][j]=hessij(p,delti,i,j);
           hess[j][i]=hess[i][j];    
   fx=func(x);          /*printf(" %lf ",hess[i][j]);*/
   for (k=1; k<=2; k++) {        }
     for (i=1;i<=npar;i++) p2[i]=x[i];      }
     p2[thetai]=x[thetai]+delti[thetai]/k;    }
     p2[thetaj]=x[thetaj]+delti[thetaj]/k;    printf("\n");
     k1=func(p2)-fx;  
       printf("\nInverting the hessian to get the covariance matrix. Wait...\n");
     p2[thetai]=x[thetai]+delti[thetai]/k;   
     p2[thetaj]=x[thetaj]-delti[thetaj]/k;    a=matrix(1,npar,1,npar);
     k2=func(p2)-fx;    y=matrix(1,npar,1,npar);
       x=vector(1,npar);
     p2[thetai]=x[thetai]-delti[thetai]/k;    indx=ivector(1,npar);
     p2[thetaj]=x[thetaj]+delti[thetaj]/k;    for (i=1;i<=npar;i++)
     k3=func(p2)-fx;      for (j=1;j<=npar;j++) a[i][j]=hess[i][j];
       ludcmp(a,npar,indx,&pd);
     p2[thetai]=x[thetai]-delti[thetai]/k;  
     p2[thetaj]=x[thetaj]-delti[thetaj]/k;    for (j=1;j<=npar;j++) {
     k4=func(p2)-fx;      for (i=1;i<=npar;i++) x[i]=0;
     res=(k1-k2-k3+k4)/4.0/delti[thetai]*k/delti[thetaj]*k/2.; /* Because of L not 2*L */      x[j]=1;
 #ifdef DEBUG      lubksb(a,npar,indx,x);
     printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);      for (i=1;i<=npar;i++){
 #endif        matcov[i][j]=x[i];
   }      }
   return res;    }
 }  
     printf("\n#Hessian matrix#\n");
 /************** Inverse of matrix **************/    for (i=1;i<=npar;i++) {
 void ludcmp(double **a, int n, int *indx, double *d)       for (j=1;j<=npar;j++) {
 {         printf("%.3e ",hess[i][j]);
   int i,imax,j,k;       }
   double big,dum,sum,temp;       printf("\n");
   double *vv;     }
    
   vv=vector(1,n);     /* Recompute Inverse */
   *d=1.0;     for (i=1;i<=npar;i++)
   for (i=1;i<=n;i++) {       for (j=1;j<=npar;j++) a[i][j]=matcov[i][j];
     big=0.0;     ludcmp(a,npar,indx,&pd);
     for (j=1;j<=n;j++)   
       if ((temp=fabs(a[i][j])) > big) big=temp;     /*  printf("\n#Hessian matrix recomputed#\n");
     if (big == 0.0) nrerror("Singular matrix in routine ludcmp");   
     vv[i]=1.0/big;     for (j=1;j<=npar;j++) {
   }       for (i=1;i<=npar;i++) x[i]=0;
   for (j=1;j<=n;j++) {       x[j]=1;
     for (i=1;i<j;i++) {       lubksb(a,npar,indx,x);
       sum=a[i][j];       for (i=1;i<=npar;i++){
       for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];         y[i][j]=x[i];
       a[i][j]=sum;         printf("%.3e ",y[i][j]);
     }       }
     big=0.0;       printf("\n");
     for (i=j;i<=n;i++) {     }
       sum=a[i][j];     */
       for (k=1;k<j;k++)   
         sum -= a[i][k]*a[k][j];     free_matrix(a,1,npar,1,npar);
       a[i][j]=sum;     free_matrix(y,1,npar,1,npar);
       if ( (dum=vv[i]*fabs(sum)) >= big) {     free_vector(x,1,npar);
         big=dum;     free_ivector(indx,1,npar);
         imax=i;     free_matrix(hess,1,npar,1,npar);
       }   
     }   
     if (j != imax) {   }
       for (k=1;k<=n;k++) {   
         dum=a[imax][k];   /*************** hessian matrix ****************/
         a[imax][k]=a[j][k];   double hessii( double x[], double delta, int theta, double delti[])
         a[j][k]=dum;   {
       }     int i;
       *d = -(*d);     int l=1, lmax=20;
       vv[imax]=vv[j];     double k1,k2;
     }     double p2[NPARMAX+1];
     indx[j]=imax;     double res;
     if (a[j][j] == 0.0) a[j][j]=TINY;     double delt, delts, nkhi=10.,nkhif=1., khi=1.e-4;
     if (j != n) {     double fx;
       dum=1.0/(a[j][j]);     int k=0,kmax=10;
       for (i=j+1;i<=n;i++) a[i][j] *= dum;     double l1;
     }   
   }     fx=func(x);
   free_vector(vv,1,n);  /* Doesn't work */    for (i=1;i<=npar;i++) p2[i]=x[i];
 ;    for(l=0 ; l <=lmax; l++){
 }       l1=pow(10,l);
       delts=delt;
 void lubksb(double **a, int n, int *indx, double b[])       for(k=1 ; k <kmax; k=k+1){
 {         delt = delta*(l1*k);
   int i,ii=0,ip,j;         p2[theta]=x[theta] +delt;
   double sum;         k1=func(p2)-fx;
          p2[theta]=x[theta]-delt;
   for (i=1;i<=n;i++) {         k2=func(p2)-fx;
     ip=indx[i];         /*res= (k1-2.0*fx+k2)/delt/delt; */
     sum=b[ip];         res= (k1+k2)/delt/delt/2.; /* Divided by because L and not 2*L */
     b[ip]=b[i];        
     if (ii)   #ifdef DEBUG
       for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];         printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx);
     else if (sum) ii=i;   #endif
     b[i]=sum;         /*if(fabs(k1-2.0*fx+k2) <1.e-13){ */
   }         if((k1 <khi/nkhi/2.) || (k2 <khi/nkhi/2.)){
   for (i=n;i>=1;i--) {           k=kmax;
     sum=b[i];         }
     for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];         else if((k1 >khi/nkhif) || (k2 >khi/nkhif)){ /* Keeps lastvalue before 3.84/2 KHI2 5% 1d.f. */
     b[i]=sum/a[i][i];           k=kmax; l=lmax*10.;
   }         }
 }         else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){
           delts=delt;
 /************ Frequencies ********************/        }
 void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx)      }
 {  /* Some frequencies */    }
      delti[theta]=delts;
   int i, m, jk;    return res;
   double ***freq; /* Frequencies */   
   double *pp;  }
   double pos;  
   FILE *ficresp;  double hessij( double x[], double delti[], int thetai,int thetaj)
   char fileresp[FILENAMELENGTH];  {
     int i;
   pp=vector(1,nlstate);    int l=1, l1, lmax=20;
     double k1,k2,k3,k4,res,fx;
   strcpy(fileresp,"p");    double p2[NPARMAX+1];
   strcat(fileresp,fileres);    int k;
   if((ficresp=fopen(fileresp,"w"))==NULL) {  
     printf("Problem with prevalence resultfile: %s\n", fileresp);    fx=func(x);
     exit(0);    for (k=1; k<=2; k++) {
   }      for (i=1;i<=npar;i++) p2[i]=x[i];
       p2[thetai]=x[thetai]+delti[thetai]/k;
   freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);      p2[thetaj]=x[thetaj]+delti[thetaj]/k;
   for (i=-1; i<=nlstate+ndeath; i++)        k1=func(p2)-fx;
     for (jk=-1; jk<=nlstate+ndeath; jk++)     
       for(m=agemin; m <= agemax+3; m++)      p2[thetai]=x[thetai]+delti[thetai]/k;
         freq[i][jk][m]=0;      p2[thetaj]=x[thetaj]-delti[thetaj]/k;
       k2=func(p2)-fx;
   for (i=1; i<=imx; i++)  {   
     for(m=firstpass; m<= lastpass-1; m++){      p2[thetai]=x[thetai]-delti[thetai]/k;
       if(agev[m][i]==0) agev[m][i]=agemax+1;      p2[thetaj]=x[thetaj]+delti[thetaj]/k;
       if(agev[m][i]==1) agev[m][i]=agemax+2;      k3=func(p2)-fx;
        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];      p2[thetai]=x[thetai]-delti[thetai]/k;
     }      p2[thetaj]=x[thetaj]-delti[thetaj]/k;
   }      k4=func(p2)-fx;
       res=(k1-k2-k3+k4)/4.0/delti[thetai]*k/delti[thetaj]*k/2.; /* Because of L not 2*L */
   fprintf(ficresp, "#");  #ifdef DEBUG
   for(i=1; i<=nlstate;i++)       printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
     fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);  #endif
 fprintf(ficresp, "\n");    }
     return res;
   for(i=(int)agemin; i <= (int)agemax+3; i++){  }
     if(i==(int)agemax+3)  
       printf("Total");  /************** Inverse of matrix **************/
     else  void ludcmp(double **a, int n, int *indx, double *d)
       printf("Age %d", i);  {
     for(jk=1; jk <=nlstate ; jk++){    int i,imax,j,k;
       for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)    double big,dum,sum,temp;
         pp[jk] += freq[jk][m][i];    double *vv;
     }   
     for(jk=1; jk <=nlstate ; jk++){    vv=vector(1,n);
       for(m=-1, pos=0; m <=0 ; m++)    *d=1.0;
         pos += freq[jk][m][i];    for (i=1;i<=n;i++) {
       if(pp[jk]>=1.e-10)      big=0.0;
         printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);      for (j=1;j<=n;j++)
       else        if ((temp=fabs(a[i][j])) > big) big=temp;
         printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);      if (big == 0.0) nrerror("Singular matrix in routine ludcmp");
     }      vv[i]=1.0/big;
     for(jk=1; jk <=nlstate ; jk++){    }
       for(m=1, pp[jk]=0; m <=nlstate+ndeath; m++)    for (j=1;j<=n;j++) {
         pp[jk] += freq[jk][m][i];      for (i=1;i<j;i++) {
     }        sum=a[i][j];
     for(jk=1,pos=0; jk <=nlstate ; jk++)        for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
       pos += pp[jk];        a[i][j]=sum;
     for(jk=1; jk <=nlstate ; jk++){      }
       if(pos>=1.e-5)      big=0.0;
         printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);      for (i=j;i<=n;i++) {
       else        sum=a[i][j];
         printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);        for (k=1;k<j;k++)
       if( i <= (int) agemax){          sum -= a[i][k]*a[k][j];
         if(pos>=1.e-5)        a[i][j]=sum;
           fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);        if ( (dum=vv[i]*fabs(sum)) >= big) {
       else          big=dum;
           fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);          imax=i;
       }        }
     }      }
     for(jk=-1; jk <=nlstate+ndeath; jk++)      if (j != imax) {
       for(m=-1; m <=nlstate+ndeath; m++)        for (k=1;k<=n;k++) {
         if(freq[jk][m][i] !=0 ) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);          dum=a[imax][k];
     if(i <= (int) agemax)          a[imax][k]=a[j][k];
       fprintf(ficresp,"\n");          a[j][k]=dum;
     printf("\n");        }
   }        *d = -(*d);
         vv[imax]=vv[j];
   fclose(ficresp);      }
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);      indx[j]=imax;
   free_vector(pp,1,nlstate);      if (a[j][j] == 0.0) a[j][j]=TINY;
       if (j != n) {
 }  /* End of Freq */        dum=1.0/(a[j][j]);
         for (i=j+1;i<=n;i++) a[i][j] *= dum;
 /************* Waves Concatenation ***************/      }
     }
 void  concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int  firstpass, int lastpass, int imx, int nlstate, int stepm)    free_vector(vv,1,n);  /* Doesn't work */
 {  ;
   /* Concatenates waves: wav[i] is the number of effective (useful waves) of individual i.  }
      Death is a valid wave (if date is known).  
      mw[mi][i] is the mi (mi=1 to wav[i])  effective wave of individual i  void lubksb(double **a, int n, int *indx, double b[])
      dh[m][i] of dh[mw[mi][i][i] is the delay between two effective waves m=mw[mi][i]  {
      and mw[mi+1][i]. dh depends on stepm.    int i,ii=0,ip,j;
      */    double sum;
    
   int i, mi, m;    for (i=1;i<=n;i++) {
   int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;      ip=indx[i];
 float sum=0.;      sum=b[ip];
       b[ip]=b[i];
   for(i=1; i<=imx; i++){      if (ii)
     mi=0;        for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
     m=firstpass;      else if (sum) ii=i;
     while(s[m][i] <= nlstate){      b[i]=sum;
       if(s[m][i]>=1)    }
         mw[++mi][i]=m;    for (i=n;i>=1;i--) {
       if(m >=lastpass)      sum=b[i];
         break;      for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];
       else      b[i]=sum/a[i][i];
         m++;    }
     }/* end while */  }
     if (s[m][i] > nlstate){  
       mi++;     /* Death is another wave */  /************ Frequencies ********************/
       /* if(mi==0)  never been interviewed correctly before death */  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)
          /* Only death is a correct wave */  {  /* Some frequencies */
       mw[mi][i]=m;   
     }    int i, m, jk, k1,i1, j1, bool, z1,z2,j;
     double ***freq; /* Frequencies */
     wav[i]=mi;    double *pp;
     if(mi==0)    double pos, k2, dateintsum=0,k2cpt=0;
       printf("Warning, no any valid information for:%d line=%d\n",num[i],i);    FILE *ficresp;
   }    char fileresp[FILENAMELENGTH];
    
   for(i=1; i<=imx; i++){    pp=vector(1,nlstate);
     for(mi=1; mi<wav[i];mi++){    probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
       if (stepm <=0)    strcpy(fileresp,"p");
         dh[mi][i]=1;    strcat(fileresp,fileres);
       else{    if((ficresp=fopen(fileresp,"w"))==NULL) {
         if (s[mw[mi+1][i]][i] > nlstate) {      printf("Problem with prevalence resultfile: %s\n", fileresp);
           j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12);       exit(0);
           if(j=0) j=1;  /* Survives at least one month after exam */    }
         }    freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
         else{    j1=0;
           j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));   
           k=k+1;    j=cptcoveff;
           if (j >= jmax) jmax=j;    if (cptcovn<1) {j=1;ncodemax[1]=1;}
           else if (j <= jmin)jmin=j;   
           sum=sum+j;    for(k1=1; k1<=j;k1++){
         }      for(i1=1; i1<=ncodemax[k1];i1++){
         jk= j/stepm;        j1++;
         jl= j -jk*stepm;        /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
         ju= j -(jk+1)*stepm;          scanf("%d", i);*/
         if(jl <= -ju)        for (i=-1; i<=nlstate+ndeath; i++)  
           dh[mi][i]=jk;          for (jk=-1; jk<=nlstate+ndeath; jk++)  
         else            for(m=agemin; m <= agemax+3; m++)
           dh[mi][i]=jk+1;              freq[i][jk][m]=0;
         if(dh[mi][i]==0)       
           dh[mi][i]=1; /* At least one step */        dateintsum=0;
       }        k2cpt=0;
     }        for (i=1; i<=imx; i++) {
   }          bool=1;
   printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,sum/k);          if  (cptcovn>0) {
 }            for (z1=1; z1<=cptcoveff; z1++)
               if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
 /*********** Health Expectancies ****************/                bool=0;
           }
 void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm)          if (bool==1) {
 {            for(m=firstpass; m<=lastpass; m++){
   /* Health expectancies */              k2=anint[m][i]+(mint[m][i]/12.);
   int i, j, nhstepm, hstepm, h;              if ((k2>=dateprev1) && (k2<=dateprev2)) {
   double age, agelim,hf;                if(agev[m][i]==0) agev[m][i]=agemax+1;
   double ***p3mat;                if(agev[m][i]==1) agev[m][i]=agemax+2;
                 if (m<lastpass) {
   FILE  *ficreseij;                  freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
   char filerese[FILENAMELENGTH];                  freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];
                 }
   strcpy(filerese,"e");               
   strcat(filerese,fileres);                if ((agev[m][i]>1) && (agev[m][i]< (agemax+3))) {
   if((ficreseij=fopen(filerese,"w"))==NULL) {                  dateintsum=dateintsum+k2;
     printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);                  k2cpt++;
   }                }
   printf("Computing Health Expectancies: result on file '%s' \n", filerese);              }
             }
   fprintf(ficreseij,"# Health expectancies\n");          }
   fprintf(ficreseij,"# Age");        }
   for(i=1; i<=nlstate;i++)         
     for(j=1; j<=nlstate;j++)        fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
       fprintf(ficreseij," %1d-%1d",i,j);  
   fprintf(ficreseij,"\n");        if  (cptcovn>0) {
           fprintf(ficresp, "\n#********** Variable ");
   hstepm=1*YEARM; /*  Every j years of age (in month) */          for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
   hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */           fprintf(ficresp, "**********\n#");
         }
   agelim=AGESUP;        for(i=1; i<=nlstate;i++)
   for (age=bage; age<=fage; age ++){ /* If stepm=6 months */          fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
     /* nhstepm age range expressed in number of stepm */        fprintf(ficresp, "\n");
     nhstepm=(int) rint((agelim-age)*YEARM/stepm);        
     /* Typically if 20 years = 20*12/6=40 stepm */         for(i=(int)agemin; i <= (int)agemax+3; i++){
     if (stepm >= YEARM) hstepm=1;          if(i==(int)agemax+3)
     nhstepm = nhstepm/hstepm;/* Expressed in hstepm, typically 40/4=10 */            printf("Total");
     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);          else
     /* Computed by stepm unit matrices, product of hstepm matrices, stored            printf("Age %d", i);
        in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */          for(jk=1; jk <=nlstate ; jk++){
     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm);              for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
               pp[jk] += freq[jk][m][i];
           }
     for(i=1; i<=nlstate;i++)          for(jk=1; jk <=nlstate ; jk++){
       for(j=1; j<=nlstate;j++)            for(m=-1, pos=0; m <=0 ; m++)
         for (h=0, eij[i][j][(int)age]=0; h<=nhstepm; h++){              pos += freq[jk][m][i];
           eij[i][j][(int)age] +=p3mat[i][j][h];            if(pp[jk]>=1.e-10)
         }              printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
                 else
     hf=1;              printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
     if (stepm >= YEARM) hf=stepm/YEARM;          }
     fprintf(ficreseij,"%.0f",age );  
     for(i=1; i<=nlstate;i++)          for(jk=1; jk <=nlstate ; jk++){
       for(j=1; j<=nlstate;j++){            for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
         fprintf(ficreseij," %.4f", hf*eij[i][j][(int)age]);              pp[jk] += freq[jk][m][i];
       }          }
     fprintf(ficreseij,"\n");  
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);          for(jk=1,pos=0; jk <=nlstate ; jk++)
   }            pos += pp[jk];
   fclose(ficreseij);          for(jk=1; jk <=nlstate ; jk++){
 }            if(pos>=1.e-5)
               printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
 /************ Variance ******************/            else
 void varevsij(char fileres[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl)              printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
 {            if( i <= (int) agemax){
   /* Variance of health expectancies */              if(pos>=1.e-5){
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/                fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);
   double **newm;                probs[i][jk][j1]= pp[jk]/pos;
   double **dnewm,**doldm;                /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/
   int i, j, nhstepm, hstepm, h;              }
   int k;              else
   FILE  *ficresvij;                fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);
   char fileresv[FILENAMELENGTH];            }
   double *xp;          }
   double **gp, **gm;         
   double ***gradg, ***trgradg;          for(jk=-1; jk <=nlstate+ndeath; jk++)
   double ***p3mat;            for(m=-1; m <=nlstate+ndeath; m++)
   double age,agelim;              if(freq[jk][m][i] !=0 ) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
   int theta;          if(i <= (int) agemax)
             fprintf(ficresp,"\n");
   strcpy(fileresv,"v");          printf("\n");
   strcat(fileresv,fileres);        }
   if((ficresvij=fopen(fileresv,"w"))==NULL) {      }
     printf("Problem with variance resultfile: %s\n", fileresv);exit(0);    }
   }    dateintmean=dateintsum/k2cpt;
   printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);   
     fclose(ficresp);
     free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
   fprintf(ficresvij,"# Covariances of life expectancies\n");    free_vector(pp,1,nlstate);
   fprintf(ficresvij,"# Age");   
   for(i=1; i<=nlstate;i++)    /* End of Freq */
     for(j=1; j<=nlstate;j++)  }
       fprintf(ficresvij," Cov(e%1d, e%1d)",i,j);  
   fprintf(ficresvij,"\n");  /************ Prevalence ********************/
   void prevalence(int agemin, float 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)
   xp=vector(1,npar);  {  /* Some frequencies */
   dnewm=matrix(1,nlstate,1,npar);   
   doldm=matrix(1,nlstate,1,nlstate);    int i, m, jk, k1, i1, j1, bool, z1,z2,j;
       double ***freq; /* Frequencies */
   hstepm=1*YEARM; /* Every year of age */    double *pp;
   hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */     double pos, k2;
   agelim = AGESUP;  
   for (age=bage; age<=fage; age ++){ /* If stepm=6 months */    pp=vector(1,nlstate);
     nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */     probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
     if (stepm >= YEARM) hstepm=1;   
     nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */    freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);    j1=0;
     gradg=ma3x(0,nhstepm,1,npar,1,nlstate);   
     gp=matrix(0,nhstepm,1,nlstate);    j=cptcoveff;
     gm=matrix(0,nhstepm,1,nlstate);    if (cptcovn<1) {j=1;ncodemax[1]=1;}
    
     for(theta=1; theta <=npar; theta++){    for(k1=1; k1<=j;k1++){
       for(i=1; i<=npar; i++){ /* Computes gradient */      for(i1=1; i1<=ncodemax[k1];i1++){
         xp[i] = x[i] + (i==theta ?delti[theta]:0);        j1++;
       }       
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm);          for (i=-1; i<=nlstate+ndeath; i++)  
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl);          for (jk=-1; jk<=nlstate+ndeath; jk++)  
       for(j=1; j<= nlstate; j++){            for(m=agemin; m <= agemax+3; m++)
         for(h=0; h<=nhstepm; h++){              freq[i][jk][m]=0;
           for(i=1, gp[h][j]=0.;i<=nlstate;i++)       
             gp[h][j] += prlim[i][i]*p3mat[i][j][h];        for (i=1; i<=imx; i++) {
         }          bool=1;
       }          if  (cptcovn>0) {
                 for (z1=1; z1<=cptcoveff; z1++)
       for(i=1; i<=npar; i++) /* Computes gradient */              if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
         xp[i] = x[i] - (i==theta ?delti[theta]:0);                bool=0;
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm);            }
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl);          if (bool==1) {
       for(j=1; j<= nlstate; j++){            for(m=firstpass; m<=lastpass; m++){
         for(h=0; h<=nhstepm; h++){              k2=anint[m][i]+(mint[m][i]/12.);
           for(i=1, gm[h][j]=0.;i<=nlstate;i++)              if ((k2>=dateprev1) && (k2<=dateprev2)) {
             gm[h][j] += prlim[i][i]*p3mat[i][j][h];                if(agev[m][i]==0) agev[m][i]=agemax+1;
         }                if(agev[m][i]==1) agev[m][i]=agemax+2;
       }                if (m<lastpass) {
       for(j=1; j<= nlstate; j++)                  if (calagedate>0)
         for(h=0; h<=nhstepm; h++){                    freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
           gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];                  else
         }                    freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
     } /* End theta */                  freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i];
                 }
     trgradg =ma3x(0,nhstepm,1,nlstate,1,npar);              }
             }
     for(h=0; h<=nhstepm; h++)          }
       for(j=1; j<=nlstate;j++)        }
         for(theta=1; theta <=npar; theta++)        for(i=(int)agemin; i <= (int)agemax+3; i++){
           trgradg[h][j][theta]=gradg[h][theta][j];          for(jk=1; jk <=nlstate ; jk++){
             for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
     for(i=1;i<=nlstate;i++)              pp[jk] += freq[jk][m][i];
       for(j=1;j<=nlstate;j++)          }
         vareij[i][j][(int)age] =0.;          for(jk=1; jk <=nlstate ; jk++){
     for(h=0;h<=nhstepm;h++){            for(m=-1, pos=0; m <=0 ; m++)
       for(k=0;k<=nhstepm;k++){              pos += freq[jk][m][i];
         matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);          }
         matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);         
         for(i=1;i<=nlstate;i++)          for(jk=1; jk <=nlstate ; jk++){
           for(j=1;j<=nlstate;j++)            for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
             vareij[i][j][(int)age] += doldm[i][j];              pp[jk] += freq[jk][m][i];
       }          }
     }         
     h=1;          for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];
     if (stepm >= YEARM) h=stepm/YEARM;         
     fprintf(ficresvij,"%.0f ",age );          for(jk=1; jk <=nlstate ; jk++){    
     for(i=1; i<=nlstate;i++)            if( i <= (int) agemax){
       for(j=1; j<=nlstate;j++){              if(pos>=1.e-5){
         fprintf(ficresvij," %.4f", h*vareij[i][j][(int)age]);                probs[i][jk][j1]= pp[jk]/pos;
       }              }
     fprintf(ficresvij,"\n");            }
     free_matrix(gp,0,nhstepm,1,nlstate);          }
     free_matrix(gm,0,nhstepm,1,nlstate);         
     free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate);        }
     free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);      }
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);    }
   } /* End age */  
   fclose(ficresvij);   
   free_vector(xp,1,npar);    free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_matrix(doldm,1,nlstate,1,npar);    free_vector(pp,1,nlstate);
   free_matrix(dnewm,1,nlstate,1,nlstate);   
   }  /* End of Freq */
 }  
   /************* Waves Concatenation ***************/
 /************ Variance of prevlim ******************/  
 void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl)  void  concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int  firstpass, int lastpass, int imx, int nlstate, int stepm)
 {  {
   /* Variance of health expectancies */    /* Concatenates waves: wav[i] is the number of effective (useful waves) of individual i.
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/       Death is a valid wave (if date is known).
   double **newm;       mw[mi][i] is the mi (mi=1 to wav[i])  effective wave of individual i
   double **dnewm,**doldm;       dh[m][i] of dh[mw[mi][i][i] is the delay between two effective waves m=mw[mi][i]
   int i, j, nhstepm, hstepm;       and mw[mi+1][i]. dh depends on stepm.
   int k;       */
   FILE  *ficresvpl;  
   char fileresvpl[FILENAMELENGTH];    int i, mi, m;
   double *xp;    /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
   double *gp, *gm;       double sum=0., jmean=0.;*/
   double **gradg, **trgradg;  
   double age,agelim;    int j, k=0,jk, ju, jl;
   int theta;    double sum=0.;
     jmin=1e+5;
   strcpy(fileresvpl,"vpl");    jmax=-1;
   strcat(fileresvpl,fileres);    jmean=0.;
   if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {    for(i=1; i<=imx; i++){
     printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);      mi=0;
     exit(0);      m=firstpass;
   }      while(s[m][i] <= nlstate){
   printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);        if(s[m][i]>=1)
           mw[++mi][i]=m;
         if(m >=lastpass)
   fprintf(ficresvpl,"# Standard deviation of prevalences limit\n");          break;
   fprintf(ficresvpl,"# Age");        else
   for(i=1; i<=nlstate;i++)          m++;
       fprintf(ficresvpl," %1d-%1d",i,i);      }/* end while */
   fprintf(ficresvpl,"\n");      if (s[m][i] > nlstate){
         mi++;     /* Death is another wave */
   xp=vector(1,npar);        /* if(mi==0)  never been interviewed correctly before death */
   dnewm=matrix(1,nlstate,1,npar);           /* Only death is a correct wave */
   doldm=matrix(1,nlstate,1,nlstate);        mw[mi][i]=m;
         }
   hstepm=1*YEARM; /* Every year of age */  
   hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */       wav[i]=mi;
   agelim = AGESUP;      if(mi==0)
   for (age=bage; age<=fage; age ++){ /* If stepm=6 months */        printf("Warning, no any valid information for:%d line=%d\n",num[i],i);
     nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */     }
     if (stepm >= YEARM) hstepm=1;  
     nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */    for(i=1; i<=imx; i++){
     gradg=matrix(1,npar,1,nlstate);      for(mi=1; mi<wav[i];mi++){
     gp=vector(1,nlstate);        if (stepm <=0)
     gm=vector(1,nlstate);          dh[mi][i]=1;
         else{
     for(theta=1; theta <=npar; theta++){          if (s[mw[mi+1][i]][i] > nlstate) {
       for(i=1; i<=npar; i++){ /* Computes gradient */            if (agedc[i] < 2*AGESUP) {
         xp[i] = x[i] + (i==theta ?delti[theta]:0);            j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12);
       }            if(j==0) j=1;  /* Survives at least one month after exam */
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl);            k=k+1;
       for(i=1;i<=nlstate;i++)            if (j >= jmax) jmax=j;
         gp[i] = prlim[i][i];            if (j <= jmin) jmin=j;
                 sum=sum+j;
       for(i=1; i<=npar; i++) /* Computes gradient */            /*if (j<0) printf("j=%d num=%d \n",j,i); */
         xp[i] = x[i] - (i==theta ?delti[theta]:0);            }
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl);          }
       for(i=1;i<=nlstate;i++)          else{
         gm[i] = prlim[i][i];            j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
             k=k+1;
       for(i=1;i<=nlstate;i++)            if (j >= jmax) jmax=j;
         gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta];            else if (j <= jmin)jmin=j;
     } /* End theta */            /*        if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
             sum=sum+j;
     trgradg =matrix(1,nlstate,1,npar);          }
           jk= j/stepm;
     for(j=1; j<=nlstate;j++)          jl= j -jk*stepm;
       for(theta=1; theta <=npar; theta++)          ju= j -(jk+1)*stepm;
         trgradg[j][theta]=gradg[theta][j];          if(jl <= -ju)
             dh[mi][i]=jk;
     for(i=1;i<=nlstate;i++)          else
       varpl[i][(int)age] =0.;            dh[mi][i]=jk+1;
     matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);          if(dh[mi][i]==0)
     matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);            dh[mi][i]=1; /* At least one step */
     for(i=1;i<=nlstate;i++)        }
       varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */      }
     }
     fprintf(ficresvpl,"%.0f ",age );    jmean=sum/k;
     for(i=1; i<=nlstate;i++)    printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
       fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age]));   }
     fprintf(ficresvpl,"\n");  /*********** Tricode ****************************/
     free_vector(gp,1,nlstate);  void tricode(int *Tvar, int **nbcode, int imx)
     free_vector(gm,1,nlstate);  {
     free_matrix(gradg,1,npar,1,nlstate);    int Ndum[20],ij=1, k, j, i;
     free_matrix(trgradg,1,nlstate,1,npar);    int cptcode=0;
   } /* End age */    cptcoveff=0;
   fclose(ficresvpl);   
   free_vector(xp,1,npar);    for (k=0; k<19; k++) Ndum[k]=0;
   free_matrix(doldm,1,nlstate,1,npar);    for (k=1; k<=7; k++) ncodemax[k]=0;
   free_matrix(dnewm,1,nlstate,1,nlstate);  
     for (j=1; j<=(cptcovn+2*cptcovprod); j++) {
 }      for (i=1; i<=imx; i++) {
         ij=(int)(covar[Tvar[j]][i]);
         Ndum[ij]++;
         /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
 /***********************************************/        if (ij > cptcode) cptcode=ij;
 /**************** Main Program *****************/      }
 /***********************************************/  
       for (i=0; i<=cptcode; i++) {
 /*int main(int argc, char *argv[])*/        if(Ndum[i]!=0) ncodemax[j]++;
 int main()      }
 {      ij=1;
   
   int i,j, k, n=MAXN,iter,m,size;  
   double agedeb, agefin,hf;      for (i=1; i<=ncodemax[j]; i++) {
   double agemin=1.e20, agemax=-1.e20;        for (k=0; k<=19; k++) {
           if (Ndum[k] != 0) {
   double fret;            nbcode[Tvar[j]][ij]=k;
   double **xi,tmp,delta;           
             ij++;
   double dum; /* Dummy variable */          }
   double ***p3mat;          if (ij > ncodemax[j]) break;
   int *indx;        }  
   char line[MAXLINE], linepar[MAXLINE];      }
   char title[MAXLINE];    }  
   char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];  
   char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];   for (k=0; k<19; k++) Ndum[k]=0;
   char filerest[FILENAMELENGTH];  
   char fileregp[FILENAMELENGTH];   for (i=1; i<=ncovmodel-2; i++) {
   char path[80],pathc[80],pathcd[80],pathtot[80];        ij=Tvar[i];
   int firstobs=1, lastobs=10;        Ndum[ij]++;
   int sdeb, sfin; /* Status at beginning and end */      }
   int c,  h , cpt,l;  
   int ju,jl, mi;   ij=1;
   int i1,j1, k1,jk,aa,bb, stepsize;   for (i=1; i<=10; i++) {
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;     if((Ndum[i]!=0) && (i<=ncovcol)){
          Tvaraff[ij]=i;
   int hstepm, nhstepm;       ij++;
   double bage, fage, age, agelim, agebase;     }
   double ftolpl=FTOL;   }
   double **prlim;   
   double *severity;      cptcoveff=ij-1;
   double ***param; /* Matrix of parameters */  }
   double  *p;  
   double **matcov; /* Matrix of covariance */  /*********** Health Expectancies ****************/
   double ***delti3; /* Scale */  
   double *delti; /* Scale */  void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm,double delti[],double **matcov )
   double ***eij, ***vareij;  
   double **varpl; /* Variances of prevalence limits by age */  {
   double *epj, vepp;    /* Health expectancies */
   char version[80]="Imach version 0.64, May 2000, INED-EUROREVES ";    int i, j, nhstepm, hstepm, h, nstepm, k, cptj;
   char *alph[]={"a","a","b","c","d","e"}, str[4];    double age, agelim, hf;
   char z[1]="c";    double ***p3mat,***varhe;
 #include <sys/time.h>    double **dnewm,**doldm;
 #include <time.h>    double *xp;
     double **gp, **gm;
   /* long total_usecs;    double ***gradg, ***trgradg;
   struct timeval start_time, end_time;    int theta;
     
   gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */    varhe=ma3x(1,nlstate*2,1,nlstate*2,(int) bage, (int) fage);
     xp=vector(1,npar);
     dnewm=matrix(1,nlstate*2,1,npar);
   printf("\nIMACH, Version 0.64");    doldm=matrix(1,nlstate*2,1,nlstate*2);
   printf("\nEnter the parameter file name: ");   
 #define windows 1    fprintf(ficreseij,"# Health expectancies\n");
 #ifdef windows    fprintf(ficreseij,"# Age");
   scanf("%s",pathtot);    for(i=1; i<=nlstate;i++)
   getcwd(pathcd, size);      for(j=1; j<=nlstate;j++)
   cut(path,optionfile,pathtot);        fprintf(ficreseij," %1d-%1d (SE)",i,j);
   chdir(path);    fprintf(ficreseij,"\n");
   replace(pathc,path);  
 #endif    if(estepm < stepm){
 #ifdef unix      printf ("Problem %d lower than %d\n",estepm, stepm);
   scanf("%s",optionfile);    }
 #endif    else  hstepm=estepm;  
     /* We compute the life expectancy from trapezoids spaced every estepm months
 /*-------- arguments in the command line --------*/     * This is mainly to measure the difference between two models: for example
      * if stepm=24 months pijx are given only every 2 years and by summing them
   strcpy(fileres,"r");     * we are calculating an estimate of the Life Expectancy assuming a linear
   strcat(fileres, optionfile);     * progression inbetween and thus overestimating or underestimating according
      * to the curvature of the survival function. If, for the same date, we
   /*---------arguments file --------*/     * estimate the model with stepm=1 month, we can keep estepm to 24 months
      * to compare the new estimate of Life expectancy with the same linear
   if((ficpar=fopen(optionfile,"r"))==NULL)    {     * hypothesis. A more precise result, taking into account a more precise
     printf("Problem with optionfile %s\n",optionfile);     * curvature will be obtained if estepm is as small as stepm. */
     goto end;  
   }    /* For example we decided to compute the life expectancy with the smallest unit */
     /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.
   strcpy(filereso,"o");       nhstepm is the number of hstepm from age to agelim
   strcat(filereso,fileres);       nstepm is the number of stepm from age to agelin.
   if((ficparo=fopen(filereso,"w"))==NULL) {       Look at hpijx to understand the reason of that which relies in memory size
     printf("Problem with Output resultfile: %s\n", filereso);goto end;       and note for a fixed period like estepm months */
   }    /* We decided (b) to get a life expectancy respecting the most precise curvature of the
        survival function given by stepm (the optimization length). Unfortunately it
 /*--------- index.htm --------*/       means that if the survival funtion is printed only each two years of age and if
        you sum them up and add 1 year (area under the trapezoids) you won't get the same
   if((fichtm=fopen("index.htm","w"))==NULL)    {       results. So we changed our mind and took the option of the best precision.
     printf("Problem with index.htm \n");goto end;    */
   }    hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */
   
  fprintf(fichtm,"<body><ul><li>Outputs files<br><br>\n    agelim=AGESUP;
         - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n    for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
 - Estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>      /* nhstepm age range expressed in number of stepm */
         - Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>      nstepm=(int) rint((agelim-age)*YEARM/stepm);
         - Transition probabilities: <a href=\"pij%s\">pij%s</a><br>      /* Typically if 20 years nstepm = 20*12/6=40 stepm */
         - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>      /* if (stepm >= YEARM) hstepm=1;*/
         - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>      nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
         - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>      p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
         - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>      gradg=ma3x(0,nhstepm,1,npar,1,nlstate*2);
         - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br><br>",fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);      gp=matrix(0,nhstepm,1,nlstate*2);
       gm=matrix(0,nhstepm,1,nlstate*2);
  fprintf(fichtm," <li>Graphs<br> <br>");  
        /* Computed by stepm unit matrices, product of hstepm matrices, stored
 for(cpt=1; cpt<nlstate;cpt++)         in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
    fprintf(fichtm,"- Prevalence of disability: p%s1.gif<br>      hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);  
 <img src=\"p%s1.gif\"><br>",strtok(optionfile, "."),strtok(optionfile, "."),cpt);   
  for(cpt=1; cpt<=nlstate;cpt++)  
      fprintf(fichtm,"- Observed and stationary  prevalence (with confident      hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
 interval) in state (%d): v%s%d.gif <br>  
 <img src=\"v%s%d.gif\"><br>",cpt,strtok(optionfile, "."),cpt,strtok(optionfile, "."),cpt);      /* Computing Variances of health expectancies */
    
  for(cpt=1; cpt<=nlstate;cpt++)       for(theta=1; theta <=npar; theta++){
      fprintf(fichtm,"- Health life expectancies by age and initial health state (%d): exp%s%d.gif <br>        for(i=1; i<=npar; i++){
 <img src=\"ex%s%d.gif\"><br>",cpt,strtok(optionfile, "."),cpt,strtok(optionfile, "."),cpt);          xp[i] = x[i] + (i==theta ?delti[theta]:0);
            }
  fprintf(fichtm,"- Total life expectancy by age and        hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
         health expectancies in states (1) and (2): e%s.gif<br>   
         <img src=\"e%s.gif\"></li> </ul></body>",strtok(optionfile, "."),strtok(optionfile, "."));        cptj=0;
         for(j=1; j<= nlstate; j++){
           for(i=1; i<=nlstate; i++){
 fclose(fichtm);            cptj=cptj+1;
             for(h=0, gp[h][cptj]=0.; h<=nhstepm-1; h++){
   /* Reads comments: lines beginning with '#' */              gp[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;
   while((c=getc(ficpar))=='#' && c!= EOF){            }
     ungetc(c,ficpar);          }
     fgets(line, MAXLINE, ficpar);        }
     puts(line);       
     fputs(line,ficparo);       
   }        for(i=1; i<=npar; i++)
   ungetc(c,ficpar);          xp[i] = x[i] - (i==theta ?delti[theta]:0);
         hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
   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\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt);       
   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\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt);        cptj=0;
   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\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt);        for(j=1; j<= nlstate; j++){
             for(i=1;i<=nlstate;i++){
   nvar=ncov-1; /* Suppressing age as a basic covariate */            cptj=cptj+1;
               for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){
   /* Read guess parameters */              gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;
   /* Reads comments: lines beginning with '#' */            }
   while((c=getc(ficpar))=='#' && c!= EOF){          }
     ungetc(c,ficpar);        }
     fgets(line, MAXLINE, ficpar);       
     puts(line);     
     fputs(line,ficparo);  
   }        for(j=1; j<= nlstate*2; j++)
   ungetc(c,ficpar);          for(h=0; h<=nhstepm-1; h++){
               gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
   param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncov);          }
     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);  /* End theta */
       printf("%1d%1d",i,j);  
       for(k=1; k<=ncov;k++){       trgradg =ma3x(0,nhstepm,1,nlstate*2,1,npar);
         fscanf(ficpar," %lf",&param[i][j][k]);  
         printf(" %lf",param[i][j][k]);       for(h=0; h<=nhstepm-1; h++)
         fprintf(ficparo," %lf",param[i][j][k]);        for(j=1; j<=nlstate*2;j++)
       }          for(theta=1; theta <=npar; theta++)
       fscanf(ficpar,"\n");          trgradg[h][j][theta]=gradg[h][theta][j];
       printf("\n");  
       fprintf(ficparo,"\n");  
     }       for(i=1;i<=nlstate*2;i++)
           for(j=1;j<=nlstate*2;j++)
   npar= (nlstate+ndeath-1)*nlstate*ncov;          varhe[i][j][(int)age] =0.;
   p=param[1][1];  
          printf("%d||",(int)age);fflush(stdout);
   /* Reads comments: lines beginning with '#' */      for(h=0;h<=nhstepm-1;h++){
   while((c=getc(ficpar))=='#' && c!= EOF){        for(k=0;k<=nhstepm-1;k++){
     ungetc(c,ficpar);          matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);
     fgets(line, MAXLINE, ficpar);          matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]);
     puts(line);          for(i=1;i<=nlstate*2;i++)
     fputs(line,ficparo);            for(j=1;j<=nlstate*2;j++)
   }              varhe[i][j][(int)age] += doldm[i][j]*hf*hf;
   ungetc(c,ficpar);        }
       }
   delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncov);  
   delti=vector(1,npar); /* Scale of each paramater (output from hesscov) */       
   for(i=1; i <=nlstate; i++){      /* Computing expectancies */
     for(j=1; j <=nlstate+ndeath-1; j++){      for(i=1; i<=nlstate;i++)
       fscanf(ficpar,"%1d%1d",&i1,&j1);        for(j=1; j<=nlstate;j++)
       printf("%1d%1d",i,j);          for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
       fprintf(ficparo,"%1d%1d",i1,j1);            eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf;
       for(k=1; k<=ncov;k++){           
         fscanf(ficpar,"%le",&delti3[i][j][k]);  /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
         printf(" %le",delti3[i][j][k]);  
         fprintf(ficparo," %le",delti3[i][j][k]);          }
       }  
       fscanf(ficpar,"\n");      fprintf(ficreseij,"%3.0f",age );
       printf("\n");      cptj=0;
       fprintf(ficparo,"\n");      for(i=1; i<=nlstate;i++)
     }        for(j=1; j<=nlstate;j++){
   }          cptj++;
   delti=delti3[1][1];          fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) );
           }
   /* Reads comments: lines beginning with '#' */      fprintf(ficreseij,"\n");
   while((c=getc(ficpar))=='#' && c!= EOF){     
     ungetc(c,ficpar);      free_matrix(gm,0,nhstepm,1,nlstate*2);
     fgets(line, MAXLINE, ficpar);      free_matrix(gp,0,nhstepm,1,nlstate*2);
     puts(line);      free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*2);
     fputs(line,ficparo);      free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar);
   }      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   ungetc(c,ficpar);    }
       free_vector(xp,1,npar);
   matcov=matrix(1,npar,1,npar);    free_matrix(dnewm,1,nlstate*2,1,npar);
   for(i=1; i <=npar; i++){    free_matrix(doldm,1,nlstate*2,1,nlstate*2);
     fscanf(ficpar,"%s",&str);    free_ma3x(varhe,1,nlstate*2,1,nlstate*2,(int) bage, (int)fage);
     printf("%s",str);  }
     fprintf(ficparo,"%s",str);  
     for(j=1; j <=i; j++){  /************ Variance ******************/
       fscanf(ficpar," %le",&matcov[i][j]);  void varevsij(char fileres[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm)
       printf(" %.5le",matcov[i][j]);  {
       fprintf(ficparo," %.5le",matcov[i][j]);    /* Variance of health expectancies */
     }    /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
     fscanf(ficpar,"\n");    double **newm;
     printf("\n");    double **dnewm,**doldm;
     fprintf(ficparo,"\n");    int i, j, nhstepm, hstepm, h, nstepm ;
   }    int k, cptcode;
   for(i=1; i <=npar; i++)    double *xp;
     for(j=i+1;j<=npar;j++)    double **gp, **gm;
       matcov[i][j]=matcov[j][i];    double ***gradg, ***trgradg;
        double ***p3mat;
   printf("\n");    double age,agelim, hf;
       int theta;
     
   if(mle==1){     fprintf(ficresvij,"# Covariances of life expectancies\n");
     /*-------- data file ----------*/    fprintf(ficresvij,"# Age");
     if((ficres =fopen(fileres,"w"))==NULL) {    for(i=1; i<=nlstate;i++)
       printf("Problem with resultfile: %s\n", fileres);goto end;      for(j=1; j<=nlstate;j++)
     }        fprintf(ficresvij," Cov(e%1d, e%1d)",i,j);
     fprintf(ficres,"#%s\n",version);    fprintf(ficresvij,"\n");
       
     if((fic=fopen(datafile,"r"))==NULL)    {    xp=vector(1,npar);
       printf("Problem with datafile: %s\n", datafile);goto end;    dnewm=matrix(1,nlstate,1,npar);
     }    doldm=matrix(1,nlstate,1,nlstate);
        
     n= lastobs;    if(estepm < stepm){
     severity = vector(1,maxwav);      printf ("Problem %d lower than %d\n",estepm, stepm);
     outcome=imatrix(1,maxwav+1,1,n);    }
     num=ivector(1,n);    else  hstepm=estepm;  
     moisnais=vector(1,n);    /* For example we decided to compute the life expectancy with the smallest unit */
     annais=vector(1,n);    /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.
     moisdc=vector(1,n);       nhstepm is the number of hstepm from age to agelim
     andc=vector(1,n);       nstepm is the number of stepm from age to agelin.
     agedc=vector(1,n);       Look at hpijx to understand the reason of that which relies in memory size
     cod=ivector(1,n);       and note for a fixed period like k years */
     weight=vector(1,n);    /* We decided (b) to get a life expectancy respecting the most precise curvature of the
     for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */       survival function given by stepm (the optimization length). Unfortunately it
     mint=matrix(1,maxwav,1,n);       means that if the survival funtion is printed only each two years of age and if
     anint=matrix(1,maxwav,1,n);       you sum them up and add 1 year (area under the trapezoids) you won't get the same
     covar=matrix(1,NCOVMAX,1,n);       results. So we changed our mind and took the option of the best precision.
     s=imatrix(1,maxwav+1,1,n);    */
     adl=imatrix(1,maxwav+1,1,n);        hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */
     tab=ivector(1,NCOVMAX);    agelim = AGESUP;
     i=1;     for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
     while (fgets(line, MAXLINE, fic) != NULL)    {      nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
       if ((i >= firstobs) && (i <lastobs)) {      nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
 sscanf(line,"%d %lf %lf %lf %lf/%lf %lf/%lf %lf/%lf %d %lf/%lf %d %lf/%lf %d %lf/%lf %d", &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]);      p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
         i=i+1;      gradg=ma3x(0,nhstepm,1,npar,1,nlstate);
       }      gp=matrix(0,nhstepm,1,nlstate);
     }       gm=matrix(0,nhstepm,1,nlstate);
   imx=i-1; /* Number of individuals */  
       for(theta=1; theta <=npar; theta++){
     fclose(fic);        for(i=1; i<=npar; i++){ /* Computes gradient */
           xp[i] = x[i] + (i==theta ?delti[theta]:0);
     if (weightopt != 1) { /* Maximisation without weights*/        }
       for(i=1;i<=n;i++) weight[i]=1.0;        hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
     }        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
     /*-calculation of age at interview from date of interview and age at death -*/  
     agev=matrix(1,maxwav,1,imx);        if (popbased==1) {
               for(i=1; i<=nlstate;i++)
     for (i=1; i<=imx; i++)  {            prlim[i][i]=probs[(int)age][i][ij];
       agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);        }
       for(m=1; (m<= maxwav); m++){   
         if(s[m][i] >0){        for(j=1; j<= nlstate; j++){
           if (s[m][i] == nlstate+1) {          for(h=0; h<=nhstepm; h++){
             if(agedc[i]>0)            for(i=1, gp[h][j]=0.;i<=nlstate;i++)
               agev[m][i]=agedc[i];              gp[h][j] += prlim[i][i]*p3mat[i][j][h];
             else{          }
               printf("Warning negative age at death: %d line:%d\n",num[i],i);        }
               agev[m][i]=-1;     
             }        for(i=1; i<=npar; i++) /* Computes gradient */
           }          xp[i] = x[i] - (i==theta ?delti[theta]:0);
           else if(s[m][i] !=9){ /* Should no more exist */        hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
             agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
             if(mint[m][i]==99 || anint[m][i]==9999)   
               agev[m][i]=1;        if (popbased==1) {
             else if(agev[m][i] <agemin){           for(i=1; i<=nlstate;i++)
               agemin=agev[m][i];            prlim[i][i]=probs[(int)age][i][ij];
               /*printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);*/        }
             }  
             else if(agev[m][i] >agemax){        for(j=1; j<= nlstate; j++){
               agemax=agev[m][i];          for(h=0; h<=nhstepm; h++){
              /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/            for(i=1, gm[h][j]=0.;i<=nlstate;i++)
             }              gm[h][j] += prlim[i][i]*p3mat[i][j][h];
             /*agev[m][i]=anint[m][i]-annais[i];*/          }
             /*   agev[m][i] = age[i]+2*m;*/        }
           }  
           else { /* =9 */        for(j=1; j<= nlstate; j++)
             agev[m][i]=1;          for(h=0; h<=nhstepm; h++){
             s[m][i]=-1;            gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
           }          }
         }      } /* End theta */
         else /*= 0 Unknown */  
           agev[m][i]=1;      trgradg =ma3x(0,nhstepm,1,nlstate,1,npar);
       }  
           for(h=0; h<=nhstepm; h++)
     }        for(j=1; j<=nlstate;j++)
     for (i=1; i<=imx; i++)  {          for(theta=1; theta <=npar; theta++)
       for(m=1; (m<= maxwav); m++){            trgradg[h][j][theta]=gradg[h][theta][j];
         if (s[m][i] > (nlstate+ndeath)) {  
           printf("Error: Wrong value in nlstate or ndeath\n");        hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
           goto end;      for(i=1;i<=nlstate;i++)
         }        for(j=1;j<=nlstate;j++)
       }          vareij[i][j][(int)age] =0.;
     }  
       for(h=0;h<=nhstepm;h++){
 printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);        for(k=0;k<=nhstepm;k++){
           matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);
     free_vector(severity,1,maxwav);          matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);
     free_imatrix(outcome,1,maxwav+1,1,n);          for(i=1;i<=nlstate;i++)
     free_vector(moisnais,1,n);            for(j=1;j<=nlstate;j++)
     free_vector(annais,1,n);              vareij[i][j][(int)age] += doldm[i][j]*hf*hf;
     free_matrix(mint,1,maxwav,1,n);        }
     free_matrix(anint,1,maxwav,1,n);      }
     free_vector(moisdc,1,n);  
     free_vector(andc,1,n);      fprintf(ficresvij,"%.0f ",age );
       for(i=1; i<=nlstate;i++)
            for(j=1; j<=nlstate;j++){
     wav=ivector(1,imx);          fprintf(ficresvij," %.4f", vareij[i][j][(int)age]);
     dh=imatrix(1,lastpass-firstpass+1,1,imx);        }
     mw=imatrix(1,lastpass-firstpass+1,1,imx);      fprintf(ficresvij,"\n");
          free_matrix(gp,0,nhstepm,1,nlstate);
     /* Concatenates waves */      free_matrix(gm,0,nhstepm,1,nlstate);
       concatwav(wav, dh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);      free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate);
           free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);
    /* Calculates basic frequencies. Computes observed prevalence at single age      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
        and prints on file fileres'p'. */    } /* End age */
       freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx);    
     free_vector(xp,1,npar);
     pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */    free_matrix(doldm,1,nlstate,1,npar);
     oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */    free_matrix(dnewm,1,nlstate,1,nlstate);
     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 */  
       /************ Variance of prevlim ******************/
     /* For Powell, parameters are in a vector p[] starting at p[1]  void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij)
        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) */    /* Variance of prevalence limit */
         /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
     mlikeli(ficres,p, npar, ncov, nlstate, ftol, func);    double **newm;
     double **dnewm,**doldm;
         int i, j, nhstepm, hstepm;
     /*--------- results files --------------*/    int k, cptcode;
     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\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, mle,weightopt);    double *xp;
         double *gp, *gm;
    jk=1;    double **gradg, **trgradg;
    fprintf(ficres,"# Parameters\n");    double age,agelim;
    printf("# Parameters\n");    int theta;
    for(i=1,jk=1; i <=nlstate; i++){     
      for(k=1; k <=(nlstate+ndeath); k++){    fprintf(ficresvpl,"# Standard deviation of prevalences limit\n");
        if (k != i)     fprintf(ficresvpl,"# Age");
          {    for(i=1; i<=nlstate;i++)
            printf("%d%d ",i,k);        fprintf(ficresvpl," %1d-%1d",i,i);
            fprintf(ficres,"%1d%1d ",i,k);    fprintf(ficresvpl,"\n");
            for(j=1; j <=ncov; j++){  
              printf("%f ",p[jk]);    xp=vector(1,npar);
              fprintf(ficres,"%f ",p[jk]);    dnewm=matrix(1,nlstate,1,npar);
              jk++;     doldm=matrix(1,nlstate,1,nlstate);
            }   
            printf("\n");    hstepm=1*YEARM; /* Every year of age */
            fprintf(ficres,"\n");    hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */
          }    agelim = AGESUP;
      }    for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
    }      nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
       if (stepm >= YEARM) hstepm=1;
     /* Computing hessian and covariance matrix */      nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
     ftolhess=ftol; /* Usually correct */      gradg=matrix(1,npar,1,nlstate);
     hesscov(matcov, p, npar, delti, ftolhess, func);      gp=vector(1,nlstate);
     fprintf(ficres,"# Scales\n");      gm=vector(1,nlstate);
     printf("# Scales\n");  
      for(i=1,jk=1; i <=nlstate; i++){      for(theta=1; theta <=npar; theta++){
       for(j=1; j <=nlstate+ndeath; j++){        for(i=1; i<=npar; i++){ /* Computes gradient */
         if (j!=i) {          xp[i] = x[i] + (i==theta ?delti[theta]:0);
           fprintf(ficres,"%1d%1d",i,j);        }
           printf("%1d%1d",i,j);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
           for(k=1; k<=ncov;k++){        for(i=1;i<=nlstate;i++)
             printf(" %.5e",delti[jk]);          gp[i] = prlim[i][i];
             fprintf(ficres," %.5e",delti[jk]);     
             jk++;        for(i=1; i<=npar; i++) /* Computes gradient */
           }          xp[i] = x[i] - (i==theta ?delti[theta]:0);
           printf("\n");        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
           fprintf(ficres,"\n");        for(i=1;i<=nlstate;i++)
         }          gm[i] = prlim[i][i];
       }  
       }        for(i=1;i<=nlstate;i++)
               gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta];
     k=1;      } /* End theta */
     fprintf(ficres,"# Covariance\n");  
     printf("# Covariance\n");      trgradg =matrix(1,nlstate,1,npar);
     for(i=1;i<=npar;i++){  
       /*  if (k>nlstate) k=1;      for(j=1; j<=nlstate;j++)
       i1=(i-1)/(ncov*nlstate)+1;         for(theta=1; theta <=npar; theta++)
       fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);          trgradg[j][theta]=gradg[theta][j];
       printf("%s%d%d",alph[k],i1,tab[i]);*/  
       fprintf(ficres,"%3d",i);      for(i=1;i<=nlstate;i++)
       printf("%3d",i);        varpl[i][(int)age] =0.;
       for(j=1; j<=i;j++){      matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);
         fprintf(ficres," %.5e",matcov[i][j]);      matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);
         printf(" %.5e",matcov[i][j]);      for(i=1;i<=nlstate;i++)
       }        varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */
       fprintf(ficres,"\n");  
       printf("\n");      fprintf(ficresvpl,"%.0f ",age );
       k++;      for(i=1; i<=nlstate;i++)
     }        fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age]));
           fprintf(ficresvpl,"\n");
     while((c=getc(ficpar))=='#' && c!= EOF){      free_vector(gp,1,nlstate);
       ungetc(c,ficpar);      free_vector(gm,1,nlstate);
       fgets(line, MAXLINE, ficpar);      free_matrix(gradg,1,npar,1,nlstate);
       puts(line);      free_matrix(trgradg,1,nlstate,1,npar);
       fputs(line,ficparo);    } /* End age */
     }  
     ungetc(c,ficpar);    free_vector(xp,1,npar);
       free_matrix(doldm,1,nlstate,1,npar);
     fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemax, &bage, &fage);    free_matrix(dnewm,1,nlstate,1,nlstate);
       
     if (fage <= 2) {  }
       bage = agemin;  
       fage = agemax;  /************ Variance of one-step probabilities  ******************/
     }  void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)
   {
     fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");    int i, j, i1, k1, j1, z1;
     fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);    int k=0, cptcode;
 /*------------ gnuplot -------------*/    double **dnewm,**doldm;
 chdir(pathcd);    double *xp;
   if((ficgp=fopen("graph.gp","w"))==NULL) {    double *gp, *gm;
     printf("Problem with file graph.gp");goto end;    double **gradg, **trgradg;
   }    double age,agelim, cov[NCOVMAX];
 #ifdef windows    int theta;
   fprintf(ficgp,"cd \"%s\" \n",pathc);    char fileresprob[FILENAMELENGTH];
 #endif  
    /* 1eme*/    strcpy(fileresprob,"prob");
     strcat(fileresprob,fileres);
   for (cpt=1; cpt<= nlstate ; cpt ++) {    if((ficresprob=fopen(fileresprob,"w"))==NULL) {
 #ifdef windows      printf("Problem with resultfile: %s\n", fileresprob);
     fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" u 1:%d \"\%%lf",agemin,fage,fileres,cpt*2);    }
 #endif    printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
 #ifdef unix   
 fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:%d \"\%%lf",agemin,fage,fileres,cpt*2);  fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n");
 #endif    fprintf(ficresprob,"# Age");
     for (i=1; i<= nlstate ; i ++) fprintf(ficgp," \%%lf (\%%lf)");    for(i=1; i<=nlstate;i++)
     fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" u 1:($%d+2*$%d) \"\%%lf",fileres,2*cpt,cpt*2+1);      for(j=1; j<=(nlstate+ndeath);j++)
     for (i=1; i<= nlstate ; i ++) fprintf(ficgp," \%%lf (\%%lf)");        fprintf(ficresprob," p%1d-%1d (SE)",i,j);
   fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"vpl%s\" u 1:($%d-2*$%d) \"\%%lf",fileres,2*cpt,2*cpt+1);   
      for (i=1; i<= nlstate ; i ++) fprintf(ficgp," \%%lf (\%%lf)");   
      fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" u 1:($%d) t\"Observed prevalence \" w l 2",fileres,2+4*(cpt-1));    fprintf(ficresprob,"\n");
 #ifdef unix  
 fprintf(ficgp,"\nset ter gif small size 400,300");  
 #endif    xp=vector(1,npar);
 fprintf(ficgp,"\nset out \"v%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt);    dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
      doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));
   }   
   /*2 eme*/    cov[1]=1;
      j=cptcoveff;
   fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",agemin,fage);    if (cptcovn<1) {j=1;ncodemax[1]=1;}
   for (i=1; i<= nlstate+1 ; i ++) {    j1=0;
 k=2*i;    for(k1=1; k1<=1;k1++){
     fprintf(ficgp,"\"t%s\" u 1:%d \"\%%lf \%%lf (\%%lf) \%%lf (\%%lf)",fileres,k);      for(i1=1; i1<=ncodemax[k1];i1++){
     for (j=1; j< nlstate ; j ++) fprintf(ficgp," \%%lf (\%%lf)");      j1++;
     if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");  
     else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);      if  (cptcovn>0) {
     fprintf(ficgp,"\"t%s\" u 1:($%d-2*$%d) \"\%%lf \%%lf (\%%lf) \%%lf (\%%lf)",fileres,k,k+1);        fprintf(ficresprob, "\n#********** Variable ");
     for (j=1; j< nlstate ; j ++) fprintf(ficgp," \%%lf (\%%lf)");        for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
     fprintf(ficgp,"\" t\"\" w l 0,");        fprintf(ficresprob, "**********\n#");
 fprintf(ficgp,"\"t%s\" u 1:($%d+2*$%d) \"\%%lf \%%lf (\%%lf) \%%lf (\%%lf)",fileres,k,k+1);      }
     for (j=1; j< nlstate ; j ++) fprintf(ficgp," \%%lf (\%%lf)");     
     if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");        for (age=bage; age<=fage; age ++){
 else fprintf(ficgp,"\" t\"\" w l 0,");          cov[2]=age;
   }           for (k=1; k<=cptcovn;k++) {
   fprintf(ficgp,"\nset out \"e%s.gif\" \nreplot\n\n",strtok(optionfile, "."));            cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];
            
   /*3eme*/          }
 for (cpt=1; cpt<= nlstate ; cpt ++) {          for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
   k=2+nlstate*(cpt-1);          for (k=1; k<=cptcovprod;k++)
     fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" u 1:%d t \"e%d1\" w l",agemin,fage,fileres,k,cpt);            cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
 for (i=1; i< nlstate ; i ++) {         
 fprintf(ficgp,",\"e%s\" u 1:%d t \"e%d%d\" w l",fileres,k+1,cpt,i+1);          gradg=matrix(1,npar,1,9);
 }           trgradg=matrix(1,9,1,npar);
 fprintf(ficgp,"\nset out \"ex%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt);          gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
 }          gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
       
 /* CV preval stat */          for(theta=1; theta <=npar; theta++){
 for (cpt=1; cpt<nlstate ; cpt ++) {            for(i=1; i<=npar; i++)
     k=3;              xp[i] = x[i] + (i==theta ?delti[theta]:0);
     fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u 2:($%d/($%d",agemin,agemax,fileres,k+cpt,k);           
     for (i=1; i< nlstate ; i ++)            pmij(pmmij,cov,ncovmodel,xp,nlstate);
       fprintf(ficgp,"+$%d",k+i);           
     fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);            k=0;
                 for(i=1; i<= (nlstate+ndeath); i++){
  l=3+(nlstate+ndeath)*cpt;              for(j=1; j<=(nlstate+ndeath);j++){
    fprintf(ficgp,",\"pij%s\" u 2:($%d/($%d",fileres,l+cpt,l);                k=k+1;
                  gp[k]=pmmij[i][j];
    for (i=1; i< nlstate ; i ++) {              }
    l=3+(nlstate+ndeath)*cpt;            }
     fprintf(ficgp,"+$%d",l+i);           
    }            for(i=1; i<=npar; i++)
   fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);              xp[i] = x[i] - (i==theta ?delti[theta]:0);
        
               pmij(pmmij,cov,ncovmodel,xp,nlstate);
   fprintf(ficgp,"set out \"p%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt);            k=0;
   }             for(i=1; i<=(nlstate+ndeath); i++){
               for(j=1; j<=(nlstate+ndeath);j++){
                 k=k+1;
   fclose(ficgp);                gm[k]=pmmij[i][j];
                  }
 chdir(path);            }
     free_matrix(agev,1,maxwav,1,imx);       
     free_ivector(wav,1,imx);            for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++)
     free_imatrix(dh,1,lastpass-firstpass+1,1,imx);              gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];  
     free_imatrix(mw,1,lastpass-firstpass+1,1,imx);          }
       
     free_imatrix(s,1,maxwav+1,1,n);          for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)
                 for(theta=1; theta <=npar; theta++)
                   trgradg[j][theta]=gradg[theta][j];
     free_ivector(num,1,n);         
     free_vector(agedc,1,n);          matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);
     free_vector(weight,1,n);          matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);
     free_matrix(covar,1,NCOVMAX,1,n);         
     fclose(ficparo);          pmij(pmmij,cov,ncovmodel,x,nlstate);
     fclose(ficres);         
   }          k=0;
           for(i=1; i<=(nlstate+ndeath); i++){
   /*________fin mle=1_________*/            for(j=1; j<=(nlstate+ndeath);j++){
                 k=k+1;
                 gm[k]=pmmij[i][j];
             }
   /* No more information from the sample is required now */          }
   /* Reads comments: lines beginning with '#' */       
   while((c=getc(ficpar))=='#' && c!= EOF){       /*printf("\n%d ",(int)age);
     ungetc(c,ficpar);       for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
     fgets(line, MAXLINE, ficpar);         printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
     puts(line);       }*/
     fputs(line,ficparo);  
   }          fprintf(ficresprob,"\n%d ",(int)age);
   ungetc(c,ficpar);  
             for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)
   fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemax, &bage, &fage);            fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i]));
   printf("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);        }
       }
   /*--------------- Prevalence limit --------------*/      free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
         free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
   strcpy(filerespl,"pl");      free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
   strcat(filerespl,fileres);      free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
   if((ficrespl=fopen(filerespl,"w"))==NULL) {    }
     printf("Problem with Prev limit resultfile: %s\n", filerespl);goto end;    free_vector(xp,1,npar);
   }    fclose(ficresprob);
   printf("Computing prevalence limit: result on file '%s' \n", filerespl);   
   fprintf(ficrespl,"#Prevalence limit\n");  }
   fprintf(ficrespl,"#Age ");  
   for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);  /******************* Printing html file ***********/
   fprintf(ficrespl,"\n");  void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \
      int lastpass, int stepm, int weightopt, char model[],\
   prlim=matrix(1,nlstate,1,nlstate);   int imx,int jmin, int jmax, double jmeanint,char optionfile[], \
   pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */   char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\
   oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */   char version[], int popforecast, int estepm ){
   newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */    int jj1, k1, i1, cpt;
   savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */    FILE *fichtm;
   oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */    /*char optionfilehtm[FILENAMELENGTH];*/
     
   agebase=agemin;    strcpy(optionfilehtm,optionfile);
   agelim=agemax;    strcat(optionfilehtm,".htm");
   ftolpl=1.e-10;    if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {
   for (age=agebase; age<=agelim; age++){      printf("Problem with %s \n",optionfilehtm), exit(0);
     prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl);    }
     fprintf(ficrespl,"%.0f",age );  
     for(i=1; i<=nlstate;i++)   fprintf(fichtm,"<body> <font size=\"2\">%s </font> <hr size=\"2\" color=\"#EC5E5E\"> \n
       fprintf(ficrespl," %.5f", prlim[i][i]);  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n
     fprintf(ficrespl,"\n");  \n
   }  Total number of observations=%d <br>\n
   fclose(ficrespl);  Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
     <hr  size=\"2\" color=\"#EC5E5E\">
   /*------------- h Pij x at various ages ------------*/   <ul><li>Outputs files<br>\n
      - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n
   strcpy(filerespij,"pij");  strcat(filerespij,fileres);   - Gnuplot file name: <a href=\"%s\">%s</a><br>\n
   if((ficrespij=fopen(filerespij,"w"))==NULL) {   - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n
     printf("Problem with Pij resultfile: %s\n", filerespij);goto end;   - Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>\n
   }   - Transition probabilities: <a href=\"pij%s\">pij%s</a><br>\n
   printf("Computing pij: result on file '%s' \n", filerespij);   - Life expectancies by age and initial health status (estepm=%2d months): <a href=\"e%s\">e%s</a> <br>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot,fileres,fileres,fileres,fileres,fileres,fileres,estepm,fileres,fileres);
   stepsize=(int) (stepm+YEARM-1)/YEARM;  
   if (stepm<=24) stepsize=2;   fprintf(fichtm,"\n
    - Parameter file with estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>\n
   agelim=AGESUP;    - Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n
   hstepm=stepsize*YEARM; /* Every year of age */   - Variances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n
   hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */    - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>\n
   for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */   - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
     nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */   
     nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */   if(popforecast==1) fprintf(fichtm,"\n
     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);   - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n
     oldm=oldms;savm=savms;   - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n
     hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm);            <br>",fileres,fileres,fileres,fileres);
     fprintf(ficrespij,"# Age");   else
     for(i=1; i<=nlstate;i++)     fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model);
       for(j=1; j<=nlstate+ndeath;j++)  fprintf(fichtm," <li>Graphs</li><p>");
         fprintf(ficrespij," %1d-%1d",i,j);  
     fprintf(ficrespij,"\n");   m=cptcoveff;
     for (h=0; h<=nhstepm; h++){   if (cptcovn < 1) {m=1;ncodemax[1]=1;}
       fprintf(ficrespij,"%.0f %.0f",agedeb, agedeb+ h*hstepm/YEARM*stepm );  
       for(i=1; i<=nlstate;i++)   jj1=0;
         for(j=1; j<=nlstate+ndeath;j++)   for(k1=1; k1<=m;k1++){
           fprintf(ficrespij," %.5f", p3mat[i][j][h]);     for(i1=1; i1<=ncodemax[k1];i1++){
       fprintf(ficrespij,"\n");         jj1++;
     }         if (cptcovn > 0) {
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);           fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
     fprintf(ficrespij,"\n");           for (cpt=1; cpt<=cptcoveff;cpt++)
   }             fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
   fclose(ficrespij);           fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
          }
   /*---------- Health expectancies and variances ------------*/         fprintf(fichtm,"<br>- Probabilities: pe%s%d.png<br>
     <img src=\"pe%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);    
   eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);         for(cpt=1; cpt<nlstate;cpt++){
   oldm=oldms;savm=savms;           fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.png<br>
   evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm);  <img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
            }
   vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);      for(cpt=1; cpt<=nlstate;cpt++) {
   oldm=oldms;savm=savms;         fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
   varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl);  interval) in state (%d): v%s%d%d.png <br>
   <img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  
   strcpy(filerest,"t");       }
   strcat(filerest,fileres);       for(cpt=1; cpt<=nlstate;cpt++) {
   if((ficrest=fopen(filerest,"w"))==NULL) {          fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br>
     printf("Problem with total LE resultfile: %s\n", filerest);goto end;  <img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
   }       }
   printf("Computing Total LEs with variances: file '%s' \n", filerest);       fprintf(fichtm,"\n<br>- Total life expectancy by age and
   fprintf(ficrest,"#Total LEs with variances: e.. (std) ");  health expectancies in states (1) and (2): e%s%d.png<br>
   for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);  <img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
   fprintf(ficrest,"\n");  fprintf(fichtm,"\n</body>");
      }
   hf=1;     }
   if (stepm >= YEARM) hf=stepm/YEARM;  fclose(fichtm);
   epj=vector(1,nlstate+1);  }
   for(age=bage; age <=fage ;age++){  
     prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl);  /******************* Gnuplot file **************/
     fprintf(ficrest," %.0f",age);  void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
     for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){  
       for(i=1, epj[j]=0.;i <=nlstate;i++) {    int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
         epj[j] += prlim[i][i]*hf*eij[i][j][(int)age];  
       }    strcpy(optionfilegnuplot,optionfilefiname);
       epj[nlstate+1] +=epj[j];    strcat(optionfilegnuplot,".gp.txt");
     }    if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
     for(i=1, vepp=0.;i <=nlstate;i++)      printf("Problem with file %s",optionfilegnuplot);
       for(j=1;j <=nlstate;j++)    }
         vepp += vareij[i][j][(int)age];  
     fprintf(ficrest," %.2f (%.2f)", epj[nlstate+1],hf*sqrt(vepp));  #ifdef windows
     for(j=1;j <=nlstate;j++){      fprintf(ficgp,"cd \"%s\" \n",pathc);
       fprintf(ficrest," %.2f (%.2f)", epj[j],hf*sqrt(vareij[j][j][(int)age]));  #endif
     }  m=pow(2,cptcoveff);
     fprintf(ficrest,"\n");   
   }   /* 1eme*/
   fclose(ficrest);    for (cpt=1; cpt<= nlstate ; cpt ++) {
   fclose(ficpar);     for (k1=1; k1<= m ; k1 ++) {
   free_vector(epj,1,nlstate+1);  
   #ifdef windows
   /*------- Variance limit prevalence------*/          fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);
   varpl=matrix(1,nlstate,(int) bage, (int) fage);  #endif
   oldm=oldms;savm=savms;  #ifdef unix
   varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl);  fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
     fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres);
     #endif
   free_matrix(varpl,1,nlstate,(int) bage, (int)fage);  
     for (i=1; i<= nlstate ; i ++) {
   free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);    if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
   free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);    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);
   free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);      for (i=1; i<= nlstate ; i ++) {
   free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);    if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
   free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);    else fprintf(ficgp," \%%*lf (\%%*lf)");
   free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);  }
       fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"vpl%s\" every :::%d::%d u 1:($2-2*$3) \"\%%lf",fileres,k1-1,k1-1);
   free_matrix(matcov,1,npar,1,npar);       for (i=1; i<= nlstate ; i ++) {
   free_vector(delti,1,npar);    if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
       else fprintf(ficgp," \%%*lf (\%%*lf)");
   free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncov);  }  
        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));
   printf("End of Imach\n");  #ifdef unix
   /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */  fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\n");
     #endif
   /* 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);*/    }
   /*------ End -----------*/    /*2 eme*/
   
  end:    for (k1=1; k1<= m ; k1 ++) {
 #ifdef windows      fprintf(ficgp,"\nset out \"e%s%d.png\" \n\n",strtok(optionfile, "."),k1);
  chdir(pathcd);      fprintf(ficgp,"set ylabel \"Years\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] ",ageminpar,fage);
 #endif      
  system("gnuplot graph.gp");      for (i=1; i<= nlstate+1 ; i ++) {
         k=2*i;
 #ifdef windows        fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:2 \"\%%lf",fileres,k1-1,k1-1);
   while (z[0] != 'q') {        for (j=1; j<= nlstate+1 ; j ++) {
     chdir(pathcd);     if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
     printf("\nType e to edit output files, c to start again, and q for exiting: ");    else fprintf(ficgp," \%%*lf (\%%*lf)");
     scanf("%s",z);  }  
     if (z[0] == 'c') system("./imach");        if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");
     else if (z[0] == 'e') {        else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);
       chdir(path);      fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",fileres,k1-1,k1-1);
       system("index.htm");        for (j=1; j<= nlstate+1 ; j ++) {
     }          if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
     else if (z[0] == 'q') exit(0);          else fprintf(ficgp," \%%*lf (\%%*lf)");
   }  }  
 #endif         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,");
       }
     }
    
     /*3eme*/
   
     for (k1=1; k1<= m ; k1 ++) {
       for (cpt=1; cpt<= nlstate ; cpt ++) {
         k=2+nlstate*(2*cpt-2);
         fprintf(ficgp,"\nset out \"exp%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
         fprintf(ficgp,"set ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt);
         /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
    for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
   fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
   fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
    for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
   fprintf(ficgp,"\" t \"e%d1\" w l",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+2*i,cpt,i+1);
   
         }
       }
     }
    
     /* CV preval stat */
       for (k1=1; k1<= m ; k1 ++) {
       for (cpt=1; cpt<nlstate ; cpt ++) {
         k=3;
         fprintf(ficgp,"set out \"p%s%d%d.png\" \n\n",strtok(optionfile, "."),cpt,k1);
         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);
   
         for (i=1; i< nlstate ; i ++)
           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);  
       }
     }  
    
     /* 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 out \"pe%s%d.png\" \n\n",strtok(optionfile, "."),jk);
        fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
        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;
            }
          }
        }
      }
      
      fclose(ficgp);
   }  /* end gnuplot */
   
   
   /*************** Moving average **************/
   void movingaverage(double agedeb, double fage,double ageminpar, double ***mobaverage){
   
     int i, cpt, cptcod;
       for (agedeb=ageminpar; agedeb<=fage; agedeb++)
         for (i=1; i<=nlstate;i++)
           for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
             mobaverage[(int)agedeb][i][cptcod]=0.;
      
       for (agedeb=ageminpar+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 ageminpar, 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(ageminpar, 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, ageminpar, 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>=(ageminpar-((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 %.f ",anproj1+cpt,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 anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,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=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM;
    
     prevalence(ageminpar, 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, ageminpar, 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,"# Age");
         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\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);  
          
           for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((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 %3.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," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]);
             }
             free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           }
         }
    
     /******/
   
         for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) {
           fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);  
           for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((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 %3.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," %15.2f", 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 ageminpar=1.e20,agemin=1.e20, agemaxpar=-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], rfileres[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,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedate;
   
     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.8d, May 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 */
     getcwd(pathcd, size);
   
     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 ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
     printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);
     fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,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);
         printf("%1d%1d",i,j);
         fprintf(ficparo,"%1d%1d",i1,j1);
         for(k=1; k<=ncovmodel;k++){
           fscanf(ficpar,"%le",&delti3[i][j][k]);
           printf(" %le",delti3[i][j][k]);
           fprintf(ficparo," %le",delti3[i][j][k]);
         }
         fscanf(ficpar,"\n");
         printf("\n");
         fprintf(ficparo,"\n");
       }
     }
     delti=delti3[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);
    
     matcov=matrix(1,npar,1,npar);
     for(i=1; i <=npar; i++){
       fscanf(ficpar,"%s",&str);
       printf("%s",str);
       fprintf(ficparo,"%s",str);
       for(j=1; j <=i; j++){
         fscanf(ficpar," %le",&matcov[i][j]);
         printf(" %.5le",matcov[i][j]);
         fprintf(ficparo," %.5le",matcov[i][j]);
       }
       fscanf(ficpar,"\n");
       printf("\n");
       fprintf(ficparo,"\n");
     }
     for(i=1; i <=npar; i++)
       for(j=i+1;j<=npar;j++)
         matcov[i][j]=matcov[j][i];
      
     printf("\n");
   
   
       /*-------- Rewriting paramater file ----------*/
        strcpy(rfileres,"r");    /* "Rparameterfile */
        strcat(rfileres,optionfilefiname);    /* Parameter file first name*/
        strcat(rfileres,".");    /* */
        strcat(rfileres,optionfilext);    /* Other files have txt extension */
       if((ficres =fopen(rfileres,"w"))==NULL) {
         printf("Problem writing new parameter file: %s\n", fileres);goto end;
       }
       fprintf(ficres,"#%s\n",version);
      
       /*-------- data file ----------*/
       if((fic=fopen(datafile,"r"))==NULL)    {
         printf("Problem with datafile: %s\n", datafile);goto end;
       }
   
       n= lastobs;
       severity = vector(1,maxwav);
       outcome=imatrix(1,maxwav+1,1,n);
       num=ivector(1,n);
       moisnais=vector(1,n);
       annais=vector(1,n);
       moisdc=vector(1,n);
       andc=vector(1,n);
       agedc=vector(1,n);
       cod=ivector(1,n);
       weight=vector(1,n);
       for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */
       mint=matrix(1,maxwav,1,n);
       anint=matrix(1,maxwav,1,n);
       s=imatrix(1,maxwav+1,1,n);
       adl=imatrix(1,maxwav+1,1,n);    
       tab=ivector(1,NCOVMAX);
       ncodemax=ivector(1,8);
   
       i=1;
       while (fgets(line, MAXLINE, fic) != NULL)    {
         if ((i >= firstobs) && (i <=lastobs)) {
          
           for (j=maxwav;j>=1;j--){
             cutv(stra, strb,line,' '); s[j][i]=atoi(strb);
             strcpy(line,stra);
             cutv(stra, strb,line,'/'); anint[j][i]=(double)(atoi(strb)); strcpy(line,stra);
             cutv(stra, strb,line,' '); mint[j][i]=(double)(atoi(strb)); strcpy(line,stra);
           }
          
           cutv(stra, strb,line,'/'); andc[i]=(double)(atoi(strb)); strcpy(line,stra);
           cutv(stra, strb,line,' '); moisdc[i]=(double)(atoi(strb)); strcpy(line,stra);
   
           cutv(stra, strb,line,'/'); annais[i]=(double)(atoi(strb)); strcpy(line,stra);
           cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra);
   
           cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra);
           for (j=ncovcol;j>=1;j--){
             cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
           }
           num[i]=atol(stra);
          
           /*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;
         }
       }
       /* printf("ii=%d", ij);
          scanf("%d",i);*/
     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 (s[4][i]==9)  s[4][i]=-1;
        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*/
     Tvar=ivector(1,15);
     Tprod=ivector(1,15);
     Tvaraff=ivector(1,15);
     Tvard=imatrix(1,15,1,2);
     Tage=ivector(1,15);      
      
     if (strlen(model) >1){
       j=0, j1=0, k1=1, k2=1;
       j=nbocc(model,'+');
       j1=nbocc(model,'*');
       cptcovn=j+1;
       cptcovprod=j1;
      
       strcpy(modelsav,model);
       if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){
         printf("Error. Non available option model=%s ",model);
         goto end;
       }
      
       for(i=(j+1); i>=1;i--){
         cutv(stra,strb,modelsav,'+');
         if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav);
         /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
         /*scanf("%d",i);*/
         if (strchr(strb,'*')) {
           cutv(strd,strc,strb,'*');
           if (strcmp(strc,"age")==0) {
             cptcovprod--;
             cutv(strb,stre,strd,'V');
             Tvar[i]=atoi(stre);
             cptcovage++;
               Tage[cptcovage]=i;
               /*printf("stre=%s ", stre);*/
           }
           else if (strcmp(strd,"age")==0) {
             cptcovprod--;
             cutv(strb,stre,strc,'V');
             Tvar[i]=atoi(stre);
             cptcovage++;
             Tage[cptcovage]=i;
           }
           else {
             cutv(strb,stre,strc,'V');
             Tvar[i]=ncovcol+k1;
             cutv(strb,strc,strd,'V');
             Tprod[k1]=i;
             Tvard[k1][1]=atoi(strc);
             Tvard[k1][2]=atoi(stre);
             Tvar[cptcovn+k2]=Tvard[k1][1];
             Tvar[cptcovn+k2+1]=Tvard[k1][2];
             for (k=1; k<=lastobs;k++)
               covar[ncovcol+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
             k1++;
             k2=k2+2;
           }
         }
         else {
           /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
          /*  scanf("%d",i);*/
         cutv(strd,strc,strb,'V');
         Tvar[i]=atoi(strc);
         }
         strcpy(modelsav,stra);  
         /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
           scanf("%d",i);*/
       }
   }
    
     /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
     printf("cptcovprod=%d ", cptcovprod);
     scanf("%d ",i);*/
       fclose(fic);
   
       /*  if(mle==1){*/
       if (weightopt != 1) { /* Maximisation without weights*/
         for(i=1;i<=n;i++) weight[i]=1.0;
       }
       /*-calculation of age at interview from date of interview and age at death -*/
       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;
          }
        if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1;
         }
       }
   
       for (i=1; i<=imx; i++)  {
         agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
         for(m=1; (m<= maxwav); m++){
           if(s[m][i] >0){
             if (s[m][i] >= nlstate+1) {
               if(agedc[i]>0)
                 if(moisdc[i]!=99 && andc[i]!=9999)
                   agev[m][i]=agedc[i];
               /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
              else {
                 if (andc[i]!=9999){
                 printf("Warning negative age at death: %d line:%d\n",num[i],i);
                 agev[m][i]=-1;
                 }
               }
             }
             else if(s[m][i] !=9){ /* Should no more exist */
               agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
               if(mint[m][i]==99 || anint[m][i]==9999)
                 agev[m][i]=1;
               else if(agev[m][i] <agemin){
                 agemin=agev[m][i];
                 /*printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);*/
               }
               else if(agev[m][i] >agemax){
                 agemax=agev[m][i];
                /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/
               }
               /*agev[m][i]=anint[m][i]-annais[i];*/
               /*   agev[m][i] = age[i]+2*m;*/
             }
             else { /* =9 */
               agev[m][i]=1;
               s[m][i]=-1;
             }
           }
           else /*= 0 Unknown */
             agev[m][i]=1;
         }
      
       }
       for (i=1; i<=imx; i++)  {
         for(m=1; (m<= maxwav); m++){
           if (s[m][i] > (nlstate+ndeath)) {
             printf("Error: Wrong value in nlstate or ndeath\n");  
             goto end;
           }
         }
       }
   
   printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
   
       free_vector(severity,1,maxwav);
       free_imatrix(outcome,1,maxwav+1,1,n);
       free_vector(moisnais,1,n);
       free_vector(annais,1,n);
       /* free_matrix(mint,1,maxwav,1,n);
          free_matrix(anint,1,maxwav,1,n);*/
       free_vector(moisdc,1,n);
       free_vector(andc,1,n);
   
      
       wav=ivector(1,imx);
       dh=imatrix(1,lastpass-firstpass+1,1,imx);
       mw=imatrix(1,lastpass-firstpass+1,1,imx);
      
       /* Concatenates waves */
         concatwav(wav, dh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
   
   
         Tcode=ivector(1,100);
         nbcode=imatrix(0,NCOVMAX,0,NCOVMAX);
         ncodemax[1]=1;
         if (cptcovn > 0) tricode(Tvar,nbcode,imx);
        
      codtab=imatrix(1,100,1,10);
      h=0;
      m=pow(2,cptcoveff);
    
      for(k=1;k<=cptcoveff; k++){
        for(i=1; i <=(m/pow(2,k));i++){
          for(j=1; j <= ncodemax[k]; j++){
            for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){
              h++;
              if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j;
              /*  printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/
            }
          }
        }
      }
      /* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]);
         codtab[1][2]=1;codtab[2][2]=2; */
      /* for(i=1; i <=m ;i++){
         for(k=1; k <=cptcovn; k++){
         printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff);
         }
         printf("\n");
         }
         scanf("%d",i);*/
      
      /* Calculates basic frequencies. Computes observed prevalence at single age
          and prints on file fileres'p'. */
   
      
      
       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 ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
    
   
      jk=1;
      fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
      printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
      for(i=1,jk=1; i <=nlstate; i++){
        for(k=1; k <=(nlstate+ndeath); k++){
          if (k != i)
            {
              printf("%d%d ",i,k);
              fprintf(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 (for hessian or gradient estimation)\n");
       printf("# Scales (for hessian or gradient estimation)\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 matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
       printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
       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);
       estepm=0;
       fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
       if (estepm==0 || estepm < stepm) estepm=stepm;
       if (fage <= 2) {
         bage = ageminpar;
         fage = agemaxpar;
       }
      
       fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
       fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
       fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
    
       while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       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);
    
   
      dateprev1=anprev1+mprev1/12.+jprev1/365.;
      dateprev2=anprev2+mprev2/12.+jprev2/365.;
   
     fscanf(ficpar,"pop_based=%d\n",&popbased);
     fprintf(ficparo,"pop_based=%d\n",popbased);  
     fprintf(ficres,"pop_based=%d\n",popbased);  
    
     while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
   
     fscanf(ficpar,"starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mov_average=%d\n",&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilav);
   fprintf(ficparo,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);
   fprintf(ficres,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);
   
   
   while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
   
     fscanf(ficpar,"popforecast=%d popfile=%s popfiledate=%lf/%lf/%lf last-popfiledate=%lf/%lf/%lf\n",&popforecast,popfile,&jpyram,&mpyram,&anpyram,&jpyram1,&mpyram1,&anpyram1);
     fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
     fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
   
    freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
   
   /*------------ gnuplot -------------*/
    printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, ageminpar,agemaxpar,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);
   
   /*--------- index.htm --------*/
   
     printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast,estepm);
   
    
     /*--------------- Prevalence limit --------------*/
    
     strcpy(filerespl,"pl");
     strcat(filerespl,fileres);
     if((ficrespl=fopen(filerespl,"w"))==NULL) {
       printf("Problem with Prev limit resultfile: %s\n", filerespl);goto end;
     }
     printf("Computing prevalence limit: result on file '%s' \n", filerespl);
     fprintf(ficrespl,"#Prevalence limit\n");
     fprintf(ficrespl,"#Age ");
     for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
     fprintf(ficrespl,"\n");
    
     prlim=matrix(1,nlstate,1,nlstate);
     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 */
     k=0;
     agebase=ageminpar;
     agelim=agemaxpar;
     ftolpl=1.e-10;
     i1=cptcoveff;
     if (cptcovn < 1){i1=1;}
   
     for(cptcov=1;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
           k=k+1;
           /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/
           fprintf(ficrespl,"\n#******");
           for(j=1;j<=cptcoveff;j++)
             fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
           fprintf(ficrespl,"******\n");
          
           for (age=agebase; age<=agelim; age++){
             prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
             fprintf(ficrespl,"%.0f",age );
             for(i=1; i<=nlstate;i++)
             fprintf(ficrespl," %.5f", prlim[i][i]);
             fprintf(ficrespl,"\n");
           }
         }
       }
     fclose(ficrespl);
   
     /*------------- h Pij x at various ages ------------*/
    
     strcpy(filerespij,"pij");  strcat(filerespij,fileres);
     if((ficrespij=fopen(filerespij,"w"))==NULL) {
       printf("Problem with Pij resultfile: %s\n", filerespij);goto end;
     }
     printf("Computing pij: result on file '%s' \n", filerespij);
    
     stepsize=(int) (stepm+YEARM-1)/YEARM;
     /*if (stepm<=24) stepsize=2;*/
   
     agelim=AGESUP;
     hstepm=stepsize*YEARM; /* Every year of age */
     hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
    
     k=0;
     for(cptcov=1;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
         k=k+1;
           fprintf(ficrespij,"\n#****** ");
           for(j=1;j<=cptcoveff;j++)
             fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
           fprintf(ficrespij,"******\n");
          
           for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
             nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
             nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
             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);  
             fprintf(ficrespij,"# Age");
             for(i=1; i<=nlstate;i++)
               for(j=1; j<=nlstate+ndeath;j++)
                 fprintf(ficrespij," %1d-%1d",i,j);
             fprintf(ficrespij,"\n");
              for (h=0; h<=nhstepm; h++){
               fprintf(ficrespij,"%d %.0f %.0f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
               for(i=1; i<=nlstate;i++)
                 for(j=1; j<=nlstate+ndeath;j++)
                   fprintf(ficrespij," %.5f", p3mat[i][j][h]);
               fprintf(ficrespij,"\n");
                }
             free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
             fprintf(ficrespij,"\n");
           }
       }
     }
   
     varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax);
   
     fclose(ficrespij);
   
   
     /*---------- Forecasting ------------------*/
     if((stepm == 1) && (strcmp(model,".")==0)){
       prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);
       if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);
     }
     else{
       erreur=108;
       printf("Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model);
     }
    
   
     /*---------- Health expectancies and variances ------------*/
   
     strcpy(filerest,"t");
     strcat(filerest,fileres);
     if((ficrest=fopen(filerest,"w"))==NULL) {
       printf("Problem with total LE resultfile: %s\n", filerest);goto end;
     }
     printf("Computing Total LEs with variances: file '%s' \n", filerest);
   
   
     strcpy(filerese,"e");
     strcat(filerese,fileres);
     if((ficreseij=fopen(filerese,"w"))==NULL) {
       printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
     }
     printf("Computing Health Expectancies: result on file '%s' \n", filerese);
   
    strcpy(fileresv,"v");
     strcat(fileresv,fileres);
     if((ficresvij=fopen(fileresv,"w"))==NULL) {
       printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
     }
     printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
     calagedate=-1;
   prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
   
     k=0;
     for(cptcov=1;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
         k=k+1;
         fprintf(ficrest,"\n#****** ");
         for(j=1;j<=cptcoveff;j++)
           fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         fprintf(ficrest,"******\n");
   
         fprintf(ficreseij,"\n#****** ");
         for(j=1;j<=cptcoveff;j++)
           fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         fprintf(ficreseij,"******\n");
   
         fprintf(ficresvij,"\n#****** ");
         for(j=1;j<=cptcoveff;j++)
           fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         fprintf(ficresvij,"******\n");
   
         eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
         oldm=oldms;savm=savms;
         evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov);  
    
         vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
         oldm=oldms;savm=savms;
          varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm);
      
   
    
         fprintf(ficrest,"#Total LEs with variances: e.. (std) ");
         for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
         fprintf(ficrest,"\n");
   
         epj=vector(1,nlstate+1);
         for(age=bage; age <=fage ;age++){
           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," %4.0f",age);
           for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
             for(i=1, epj[j]=0.;i <=nlstate;i++) {
               epj[j] += prlim[i][i]*eij[i][j][(int)age];
               /*  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
             }
             epj[nlstate+1] +=epj[j];
           }
   
           for(i=1, vepp=0.;i <=nlstate;i++)
             for(j=1;j <=nlstate;j++)
               vepp += vareij[i][j][(int)age];
           fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
           for(j=1;j <=nlstate;j++){
             fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
           }
           fprintf(ficrest,"\n");
         }
       }
     }
   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);
     fclose(ficreseij);
     fclose(ficresvij);
     fclose(ficrest);
     fclose(ficpar);
     free_vector(epj,1,nlstate+1);
    
     /*------- Variance limit prevalence------*/  
   
     strcpy(fileresvpl,"vpl");
     strcat(fileresvpl,fileres);
     if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
       printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);
       exit(0);
     }
     printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);
   
     k=0;
     for(cptcov=1;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
         k=k+1;
         fprintf(ficresvpl,"\n#****** ");
         for(j=1;j<=cptcoveff;j++)
           fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         fprintf(ficresvpl,"******\n");
        
         varpl=matrix(1,nlstate,(int) bage, (int) fage);
         oldm=oldms;savm=savms;
        varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);
       }
    }
   
     fclose(ficresvpl);
   
     /*---------- End : free ----------------*/
     free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
    
     free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
     free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
    
    
     free_matrix(pmmij,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(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
    
     free_matrix(matcov,1,npar,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);
   
     if(erreur >0)
       printf("End of Imach with error or warning %d\n",erreur);
     else   printf("End of Imach\n");
     /*  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 uSec.\n", total_usecs);*/
     /*------ End -----------*/
   
   
    end:
   #ifdef windows
     /* chdir(pathcd);*/
   #endif
    /*system("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
     while (z[0] != 'q') {
       /* chdir(path); */
       printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");
       scanf("%s",z);
       if (z[0] == 'c') system("./imach");
       else if (z[0] == 'e') system(optionfilehtm);
       else if (z[0] == 'g') system(plotcmd);
       else if (z[0] == 'q') exit(0);
     }
   #endif
   }
   
   

Removed from v.1.1.1.1  
changed lines
  Added in v.1.42


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