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

Diff for /imach/src/imach.c between versions 1.17 and 1.41.2.1

version 1.17, 2002/02/20 17:15:02 version 1.41.2.1, 2003/06/12 10:43:20
Line 1 Line 1
      /* $Id$
 /*********************** Imach **************************************             Interpolated Markov Chain
   This program computes Healthy Life Expectancies from cross-longitudinal  
   data. Cross-longitudinal consist in a first survey ("cross") where    Short summary of the programme:
   individuals from different ages are interviewed on their health status   
   or degree of  disability. At least a second wave of interviews    This program computes Healthy Life Expectancies from
   ("longitudinal") should  measure each new individual health status.    cross-longitudinal data. Cross-longitudinal data consist in: -1- a
   Health expectancies are computed from the transistions observed between    first survey ("cross") where individuals from different ages are
   waves and are computed for each degree of severity of disability (number    interviewed on their health status or degree of disability (in the
   of life states). More degrees you consider, more time is necessary to    case of a health survey which is our main interest) -2- at least a
   reach the Maximum Likelihood of the parameters involved in the model.    second wave of interviews ("longitudinal") which measure each change
   The simplest model is the multinomial logistic model where pij is    (if any) in individual health status.  Health expectancies are
   the probabibility to be observed in state j at the second wave conditional    computed from the time spent in each health state according to a
   to be observed in state i at the first wave. Therefore the model is:    model. More health states you consider, more time is necessary to reach the
   log(pij/pii)= aij + bij*age+ cij*sex + etc , where 'age' is age and 'sex'    Maximum Likelihood of the parameters involved in the model.  The
   is a covariate. If you want to have a more complex model than "constant and    simplest model is the multinomial logistic model where pij is the
   age", you should modify the program where the markup    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
     conditional to the observed state i at age x. The delay 'h' can be
     split into an exact number (nh*stepm) of unobserved intermediate
     states. This elementary transition (by month or quarter trimester,
     semester or year) is model as a multinomial logistic.  The hPx
     matrix is simply the matrix product of nh*stepm elementary matrices
     and the contribution of each individual to the likelihood is simply
     hPijx.
   
   Also this programme outputs the covariance matrix of the parameters but also    Also this programme outputs the covariance matrix of the parameters but also
   of the life expectancies. It also computes the prevalence limits.    of the life expectancies. It also computes the prevalence limits.
Line 48 Line 56
 #include <unistd.h>  #include <unistd.h>
   
 #define MAXLINE 256  #define MAXLINE 256
   #define GNUPLOTPROGRAM "wgnuplot"
   /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
 #define FILENAMELENGTH 80  #define FILENAMELENGTH 80
 /*#define DEBUG*/  /*#define DEBUG*/
 #define windows  
   /*#define windows*/
 #define GLOCK_ERROR_NOPATH              -1      /* empty path */  #define GLOCK_ERROR_NOPATH              -1      /* empty path */
 #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */  #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */
   
Line 67 Line 78
 #define AGEBASE 40  #define AGEBASE 40
   
   
   int erreur; /* Error number */
 int nvar;  int nvar;
 int cptcovn, cptcovage=0, cptcoveff=0,cptcov;  int cptcovn, cptcovage=0, cptcoveff=0,cptcov;
 int npar=NPARMAX;  int npar=NPARMAX;
 int nlstate=2; /* Number of live states */  int nlstate=2; /* Number of live states */
 int ndeath=1; /* Number of dead states */  int ndeath=1; /* Number of dead states */
 int ncovmodel, ncov;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */  int ncovmodel, ncovcol;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
 int popbased=0;  int popbased=0;
   
 int *wav; /* Number of waves for this individuual 0 is possible */  int *wav; /* Number of waves for this individuual 0 is possible */
Line 84  int **dh; /* dh[mi][i] is number of step Line 96  int **dh; /* dh[mi][i] is number of step
 double jmean; /* Mean space between 2 waves */  double jmean; /* Mean space between 2 waves */
 double **oldm, **newm, **savm; /* Working pointers to matrices */  double **oldm, **newm, **savm; /* Working pointers to matrices */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */  double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
 FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf;  FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
 FILE *ficgp, *fichtm,*ficresprob,*ficpop;  FILE *ficgp,*ficresprob,*ficpop;
 FILE *ficreseij;  FILE *ficreseij;
   char filerese[FILENAMELENGTH];    char filerese[FILENAMELENGTH];
  FILE  *ficresvij;   FILE  *ficresvij;
Line 113  FILE *ficreseij; Line 125  FILE *ficreseij;
 static double maxarg1,maxarg2;  static double maxarg1,maxarg2;
 #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)>(maxarg2)? (maxarg1):(maxarg2))  #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)>(maxarg2)? (maxarg1):(maxarg2))
 #define FMIN(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)<(maxarg2)? (maxarg1):(maxarg2))  #define FMIN(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1)<(maxarg2)? (maxarg1):(maxarg2))
     
 #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a))  #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a))
 #define rint(a) floor(a+0.5)  #define rint(a) floor(a+0.5)
   
Line 125  int imx; Line 137  int imx;
 int stepm;  int stepm;
 /* Stepm, step in month: minimum step interpolation*/  /* Stepm, step in month: minimum step interpolation*/
   
   int estepm;
   /* Estepm, step in month to interpolate survival function in order to approximate Life Expectancy*/
   
 int m,nb;  int m,nb;
 int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;  int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;
 double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;  double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
 double **pmmij, ***probs, ***mobaverage;  double **pmmij, ***probs, ***mobaverage;
   double dateintmean=0;
   
 double *weight;  double *weight;
 int **s; /* Status */  int **s; /* Status */
Line 139  double ftol=FTOL; /* Tolerance for compu Line 155  double ftol=FTOL; /* Tolerance for compu
 double ftolhess; /* Tolerance for computing hessian */  double ftolhess; /* Tolerance for computing hessian */
   
 /**************** split *************************/  /**************** split *************************/
 static  int split( char *path, char *dirc, char *name )  static  int split( char *path, char *dirc, char *name, char *ext, char *finame )
 {  {
    char *s;                             /* pointer */     char *s;                             /* pointer */
    int  l1, l2;                         /* length counters */     int  l1, l2;                         /* length counters */
   
    l1 = strlen( path );                 /* length of path */     l1 = strlen( path );                 /* length of path */
    if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );     if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
   #ifdef windows
    s = strrchr( path, '\\' );           /* find last / */     s = strrchr( path, '\\' );           /* find last / */
   #else
      s = strrchr( path, '/' );            /* find last / */
   #endif
    if ( s == NULL ) {                   /* no directory, so use current */     if ( s == NULL ) {                   /* no directory, so use current */
 #if     defined(__bsd__)                /* get current working directory */  #if     defined(__bsd__)                /* get current working directory */
       extern char       *getwd( );        extern char       *getwd( );
Line 169  static int split( char *path, char *dirc Line 189  static int split( char *path, char *dirc
       dirc[l1-l2] = 0;                  /* add zero */        dirc[l1-l2] = 0;                  /* add zero */
    }     }
    l1 = strlen( dirc );                 /* length of directory */     l1 = strlen( dirc );                 /* length of directory */
   #ifdef windows
    if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }     if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
   #else
      if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }
   #endif
      s = strrchr( name, '.' );            /* find last / */
      s++;
      strcpy(ext,s);                       /* save extension */
      l1= strlen( name);
      l2= strlen( s)+1;
      strncpy( finame, name, l1-l2);
      finame[l1-l2]= 0;
    return( 0 );                         /* we're done */     return( 0 );                         /* we're done */
 }  }
   
Line 660  double **prevalim(double **prlim, int nl Line 691  double **prevalim(double **prlim, int nl
     
       for (k=1; k<=cptcovn;k++) {        for (k=1; k<=cptcovn;k++) {
         cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];          cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
         /*printf("ij=%d Tvar[k]=%d nbcode=%d cov=%lf\n",ij, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k]);*/          /*      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]]);*/
       }        }
       for (k=1; k<=cptcovage;k++)        for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
         cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];  
       for (k=1; k<=cptcovprod;k++)        for (k=1; k<=cptcovprod;k++)
         cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];          cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
   
       /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/        /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
       /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/        /*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]);*/
     out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);      out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
   
     savm=oldm;      savm=oldm;
Line 717  double **pmij(double **ps, double *cov, Line 747  double **pmij(double **ps, double *cov,
         s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];          s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
         /*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/          /*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
       }        }
       ps[i][j]=(s2);        ps[i][j]=s2;
     }      }
   }    }
     /*ps[3][2]=1;*/      /*ps[3][2]=1;*/
Line 900  void mlikeli(FILE *ficres,double p[], in Line 930  void mlikeli(FILE *ficres,double p[], in
   powell(p,xi,npar,ftol,&iter,&fret,func);    powell(p,xi,npar,ftol,&iter,&fret,func);
   
    printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));     printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
   fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f ",iter,func(p));    fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
   
 }  }
   
Line 1150  void lubksb(double **a, int n, int *indx Line 1180  void lubksb(double **a, int n, int *indx
 }  }
   
 /************ Frequencies ********************/  /************ Frequencies ********************/
 void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax, int fprev1,int lprev1,double **mint,double **anint)  void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2)
 {  /* Some frequencies */  {  /* Some frequencies */
     
   int i, m, jk, k1, k2,i1, j1, bool, z1,z2,j;    int i, m, jk, k1,i1, j1, bool, z1,z2,j;
   double ***freq; /* Frequencies */    double ***freq; /* Frequencies */
   double *pp;    double *pp;
   double pos;    double pos, k2, dateintsum=0,k2cpt=0;
   FILE *ficresp;    FILE *ficresp;
   char fileresp[FILENAMELENGTH];    char fileresp[FILENAMELENGTH];
    
   pp=vector(1,nlstate);    pp=vector(1,nlstate);
   probs= ma3x(1,130 ,1,8, 1,8);    probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
   strcpy(fileresp,"p");    strcpy(fileresp,"p");
   strcat(fileresp,fileres);    strcat(fileresp,fileres);
   if((ficresp=fopen(fileresp,"w"))==NULL) {    if((ficresp=fopen(fileresp,"w"))==NULL) {
Line 1170  void  freqsummary(char fileres[], int ag Line 1200  void  freqsummary(char fileres[], int ag
   }    }
   freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);    freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
   j1=0;    j1=0;
    
   j=cptcoveff;    j=cptcoveff;
   if (cptcovn<1) {j=1;ncodemax[1]=1;}    if (cptcovn<1) {j=1;ncodemax[1]=1;}
    
   for(k1=1; k1<=j;k1++){    for(k1=1; k1<=j;k1++){
    for(i1=1; i1<=ncodemax[k1];i1++){      for(i1=1; i1<=ncodemax[k1];i1++){
        j1++;        j1++;
        /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);        /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
          scanf("%d", i);*/          scanf("%d", i);*/
         for (i=-1; i<=nlstate+ndeath; i++)          for (i=-1; i<=nlstate+ndeath; i++)  
          for (jk=-1; jk<=nlstate+ndeath; jk++)            for (jk=-1; jk<=nlstate+ndeath; jk++)  
            for(m=agemin; m <= agemax+3; m++)            for(m=agemin; m <= agemax+3; m++)
              freq[i][jk][m]=0;              freq[i][jk][m]=0;
               
        for (i=1; i<=imx; i++) {        dateintsum=0;
          bool=1;        k2cpt=0;
          if  (cptcovn>0) {        for (i=1; i<=imx; i++) {
            for (z1=1; z1<=cptcoveff; z1++)          bool=1;
              if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])          if  (cptcovn>0) {
                bool=0;            for (z1=1; z1<=cptcoveff; z1++)
          }              if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
           if (bool==1) {                bool=0;
             for(m=fprev1; m<=lprev1; m++){          }
              k2=anint[m][i]+(mint[m][i]/12.);          if (bool==1) {
              if ((k2>=1984) && (k2<=1988.5)) {            for(m=firstpass; m<=lastpass; m++){
              if(agev[m][i]==0) agev[m][i]=agemax+1;              k2=anint[m][i]+(mint[m][i]/12.);
              if(agev[m][i]==1) agev[m][i]=agemax+2;              if ((k2>=dateprev1) && (k2<=dateprev2)) {
              freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];                if(agev[m][i]==0) agev[m][i]=agemax+1;
              freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];                if(agev[m][i]==1) agev[m][i]=agemax+2;
              }                if (m<lastpass) {
                   freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
                   freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];
                 }
                
                 if ((agev[m][i]>1) && (agev[m][i]< (agemax+3))) {
                   dateintsum=dateintsum+k2;
                   k2cpt++;
                 }
             }              }
           }            }
        }          }
         if  (cptcovn>0) {        }
          fprintf(ficresp, "\n#********** Variable ");  
          for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);  
        fprintf(ficresp, "**********\n#");  
         }  
        for(i=1; i<=nlstate;i++)  
          fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);  
        fprintf(ficresp, "\n");  
                 
   for(i=(int)agemin; i <= (int)agemax+3; i++){        fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
     if(i==(int)agemax+3)  
       printf("Total");  
     else  
       printf("Age %d", i);  
     for(jk=1; jk <=nlstate ; jk++){  
       for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)  
         pp[jk] += freq[jk][m][i];  
     }  
     for(jk=1; jk <=nlstate ; jk++){  
       for(m=-1, pos=0; m <=0 ; m++)  
         pos += freq[jk][m][i];  
       if(pp[jk]>=1.e-10)  
         printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);  
       else  
         printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);  
     }  
   
      for(jk=1; jk <=nlstate ; jk++){        if  (cptcovn>0) {
       for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)          fprintf(ficresp, "\n#********** Variable ");
         pp[jk] += freq[jk][m][i];          for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
      }          fprintf(ficresp, "**********\n#");
         }
         for(i=1; i<=nlstate;i++)
           fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
         fprintf(ficresp, "\n");
        
         for(i=(int)agemin; i <= (int)agemax+3; i++){
           if(i==(int)agemax+3)
             printf("Total");
           else
             printf("Age %d", i);
           for(jk=1; jk <=nlstate ; jk++){
             for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
               pp[jk] += freq[jk][m][i];
           }
           for(jk=1; jk <=nlstate ; jk++){
             for(m=-1, pos=0; m <=0 ; m++)
               pos += freq[jk][m][i];
             if(pp[jk]>=1.e-10)
               printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
             else
               printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
           }
   
     for(jk=1,pos=0; jk <=nlstate ; jk++)          for(jk=1; jk <=nlstate ; jk++){
       pos += pp[jk];            for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
     for(jk=1; jk <=nlstate ; jk++){              pp[jk] += freq[jk][m][i];
       if(pos>=1.e-5)  
         printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);  
       else  
         printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);  
       if( i <= (int) agemax){  
         if(pos>=1.e-5){  
           fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);  
           probs[i][jk][j1]= pp[jk]/pos;  
           /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/  
         }          }
       else  
           fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);          for(jk=1,pos=0; jk <=nlstate ; jk++)
             pos += pp[jk];
           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);
             else
               printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
             if( i <= (int) agemax){
               if(pos>=1.e-5){
                 fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);
                 probs[i][jk][j1]= pp[jk]/pos;
                 /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/
               }
               else
                 fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);
             }
           }
          
           for(jk=-1; jk <=nlstate+ndeath; jk++)
             for(m=-1; m <=nlstate+ndeath; m++)
               if(freq[jk][m][i] !=0 ) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
           if(i <= (int) agemax)
             fprintf(ficresp,"\n");
           printf("\n");
       }        }
     }      }
     for(jk=-1; jk <=nlstate+ndeath; jk++)    }
       for(m=-1; m <=nlstate+ndeath; m++)    dateintmean=dateintsum/k2cpt;
         if(freq[jk][m][i] !=0 ) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);  
     if(i <= (int) agemax)  
       fprintf(ficresp,"\n");  
     printf("\n");  
     }  
     }  
  }  
     
   fclose(ficresp);    fclose(ficresp);
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);    free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);    free_vector(pp,1,nlstate);
    
 }  /* End of Freq */    /* End of Freq */
   }
   
 /************ Prevalence ********************/  /************ Prevalence ********************/
 void prevalence(int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax, int fprev1,int lprev1)  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)
 {  /* Some frequencies */  {  /* Some frequencies */
     
   int i, m, jk, k1, i1, j1, bool, z1,z2,j;    int i, m, jk, k1, i1, j1, bool, z1,z2,j;
   double ***freq; /* Frequencies */    double ***freq; /* Frequencies */
   double *pp;    double *pp;
   double pos;    double pos, k2;
   
   pp=vector(1,nlstate);    pp=vector(1,nlstate);
   probs= ma3x(1,130 ,1,8, 1,8);    probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
     
   freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);    freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
   j1=0;    j1=0;
Line 1293  void prevalence(int agemin, int agemax, Line 1338  void prevalence(int agemin, int agemax,
       for (i=-1; i<=nlstate+ndeath; i++)          for (i=-1; i<=nlstate+ndeath; i++)  
         for (jk=-1; jk<=nlstate+ndeath; jk++)            for (jk=-1; jk<=nlstate+ndeath; jk++)  
           for(m=agemin; m <= agemax+3; m++)            for(m=agemin; m <= agemax+3; m++)
           freq[i][jk][m]=0;              freq[i][jk][m]=0;
             
       for (i=1; i<=imx; i++) {        for (i=1; i<=imx; i++) {
         bool=1;          bool=1;
         if  (cptcovn>0) {          if  (cptcovn>0) {
           for (z1=1; z1<=cptcoveff; z1++)            for (z1=1; z1<=cptcoveff; z1++)
             if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])              if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
               bool=0;                bool=0;
               }          }
         if (bool==1) {          if (bool==1) {
           for(m=fprev1; m<=lprev1; m++){            for(m=firstpass; m<=lastpass; m++){
             if(agev[m][i]==0) agev[m][i]=agemax+1;              k2=anint[m][i]+(mint[m][i]/12.);
             if(agev[m][i]==1) agev[m][i]=agemax+2;              if ((k2>=dateprev1) && (k2<=dateprev2)) {
             freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];                if(agev[m][i]==0) agev[m][i]=agemax+1;
             freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];                if(agev[m][i]==1) agev[m][i]=agemax+2;
                 if (m<lastpass)
                   if (calagedate>0) freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
                 else
                  freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
                  freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i];
               }
           }            }
         }          }
       }        }
        for(i=(int)agemin; i <= (int)agemax+3; i++){          for(i=(int)agemin; i <= (int)agemax+3; i++){
         for(jk=1; jk <=nlstate ; jk++){            for(jk=1; jk <=nlstate ; jk++){
           for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)              for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
             pp[jk] += freq[jk][m][i];                pp[jk] += freq[jk][m][i];
         }            }
         for(jk=1; jk <=nlstate ; jk++){            for(jk=1; jk <=nlstate ; jk++){
           for(m=-1, pos=0; m <=0 ; m++)              for(m=-1, pos=0; m <=0 ; m++)
             pos += freq[jk][m][i];              pos += freq[jk][m][i];
         }          }
                 
Line 1336  void prevalence(int agemin, int agemax, Line 1387  void prevalence(int agemin, int agemax,
            }             }
          }           }
                     
          }          }
     }      }
   }    }
    
     
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);    free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);    free_vector(pp,1,nlstate);
     
 }  /* End of Freq */  }  /* End of Freq */
   
 /************* Waves Concatenation ***************/  /************* 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)  void  concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int  firstpass, int lastpass, int imx, int nlstate, int stepm)
Line 1401  void  concatwav(int wav[], int **dh, int Line 1453  void  concatwav(int wav[], int **dh, int
           if (j >= jmax) jmax=j;            if (j >= jmax) jmax=j;
           if (j <= jmin) jmin=j;            if (j <= jmin) jmin=j;
           sum=sum+j;            sum=sum+j;
           /* if (j<10) printf("j=%d num=%d ",j,i); */            /*if (j<0) printf("j=%d num=%d \n",j,i); */
           }            }
         }          }
         else{          else{
Line 1409  void  concatwav(int wav[], int **dh, int Line 1461  void  concatwav(int wav[], int **dh, int
           k=k+1;            k=k+1;
           if (j >= jmax) jmax=j;            if (j >= jmax) jmax=j;
           else if (j <= jmin)jmin=j;            else if (j <= jmin)jmin=j;
           /*   if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */            /*        if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
           sum=sum+j;            sum=sum+j;
         }          }
         jk= j/stepm;          jk= j/stepm;
Line 1455  void tricode(int *Tvar, int **nbcode, in Line 1507  void tricode(int *Tvar, int **nbcode, in
       for (k=0; k<=19; k++) {        for (k=0; k<=19; k++) {
         if (Ndum[k] != 0) {          if (Ndum[k] != 0) {
           nbcode[Tvar[j]][ij]=k;            nbcode[Tvar[j]][ij]=k;
            
           ij++;            ij++;
         }          }
         if (ij > ncodemax[j]) break;          if (ij > ncodemax[j]) break;
Line 1471  void tricode(int *Tvar, int **nbcode, in Line 1524  void tricode(int *Tvar, int **nbcode, in
   
  ij=1;   ij=1;
  for (i=1; i<=10; i++) {   for (i=1; i<=10; i++) {
    if((Ndum[i]!=0) && (i<=ncov)){     if((Ndum[i]!=0) && (i<=ncovcol)){
      Tvaraff[ij]=i;       Tvaraff[ij]=i;
      ij++;       ij++;
    }     }
Line 1482  void tricode(int *Tvar, int **nbcode, in Line 1535  void tricode(int *Tvar, int **nbcode, in
   
 /*********** Health Expectancies ****************/  /*********** Health Expectancies ****************/
   
 void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij)  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 )
   
 {  {
   /* Health expectancies */    /* Health expectancies */
   int i, j, nhstepm, hstepm, h;    int i, j, nhstepm, hstepm, h, nstepm, k, cptj;
   double age, agelim,hf;    double age, agelim, hf;
   double ***p3mat;    double ***p3mat,***varhe;
     double **dnewm,**doldm;
     double *xp;
     double **gp, **gm;
     double ***gradg, ***trgradg;
     int theta;
   
     varhe=ma3x(1,nlstate*2,1,nlstate*2,(int) bage, (int) fage);
     xp=vector(1,npar);
     dnewm=matrix(1,nlstate*2,1,npar);
     doldm=matrix(1,nlstate*2,1,nlstate*2);
     
   fprintf(ficreseij,"# Health expectancies\n");    fprintf(ficreseij,"# Health expectancies\n");
   fprintf(ficreseij,"# Age");    fprintf(ficreseij,"# Age");
   for(i=1; i<=nlstate;i++)    for(i=1; i<=nlstate;i++)
     for(j=1; j<=nlstate;j++)      for(j=1; j<=nlstate;j++)
       fprintf(ficreseij," %1d-%1d",i,j);        fprintf(ficreseij," %1d-%1d (SE)",i,j);
   fprintf(ficreseij,"\n");    fprintf(ficreseij,"\n");
   
   hstepm=1*YEARM; /*  Every j years of age (in month) */    if(estepm < stepm){
   hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */      printf ("Problem %d lower than %d\n",estepm, stepm);
     }
     else  hstepm=estepm;  
     /* We compute the life expectancy from trapezoids spaced every estepm months
      * 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
      * we are calculating an estimate of the Life Expectancy assuming a linear
      * progression inbetween and thus overestimating or underestimating according
      * to the curvature of the survival function. If, for the same date, we
      * 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
      * hypothesis. A more precise result, taking into account a more precise
      * curvature will be obtained if estepm is as small as stepm. */
   
     /* 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.
        nhstepm is the number of hstepm from age to agelim
        nstepm is the number of stepm from age to agelin.
        Look at hpijx to understand the reason of that which relies in memory size
        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
        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
        results. So we changed our mind and took the option of the best precision.
     */
     hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */
   
   agelim=AGESUP;    agelim=AGESUP;
   for (age=bage; age<=fage; age ++){ /* If stepm=6 months */    for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
     /* nhstepm age range expressed in number of stepm */      /* nhstepm age range expressed in number of stepm */
     nhstepm=(int) rint((agelim-age)*YEARM/stepm);      nstepm=(int) rint((agelim-age)*YEARM/stepm);
     /* Typically if 20 years = 20*12/6=40 stepm */      /* Typically if 20 years nstepm = 20*12/6=40 stepm */
     if (stepm >= YEARM) hstepm=1;      /* if (stepm >= YEARM) hstepm=1;*/
     nhstepm = nhstepm/hstepm;/* Expressed in hstepm, typically 40/4=10 */      nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
       gradg=ma3x(0,nhstepm,1,npar,1,nlstate*2);
       gp=matrix(0,nhstepm,1,nlstate*2);
       gm=matrix(0,nhstepm,1,nlstate*2);
   
     /* Computed by stepm unit matrices, product of hstepm matrices, stored      /* Computed by stepm unit matrices, product of hstepm matrices, stored
        in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */         in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);        hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);  
    
   
       hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
   
       /* Computing Variances of health expectancies */
   
        for(theta=1; theta <=npar; theta++){
         for(i=1; i<=npar; i++){
           xp[i] = x[i] + (i==theta ?delti[theta]:0);
         }
         hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
    
         cptj=0;
         for(j=1; j<= nlstate; j++){
           for(i=1; i<=nlstate; i++){
             cptj=cptj+1;
             for(h=0, gp[h][cptj]=0.; h<=nhstepm-1; h++){
               gp[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;
             }
           }
         }
        
        
         for(i=1; i<=npar; i++)
           xp[i] = x[i] - (i==theta ?delti[theta]:0);
         hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
        
         cptj=0;
         for(j=1; j<= nlstate; j++){
           for(i=1;i<=nlstate;i++){
             cptj=cptj+1;
             for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){
               gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;
             }
           }
         }
        
      
   
         for(j=1; j<= nlstate*2; j++)
           for(h=0; h<=nhstepm-1; h++){
             gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
           }
   
        }
      
   /* End theta */
   
        trgradg =ma3x(0,nhstepm,1,nlstate*2,1,npar);
   
        for(h=0; h<=nhstepm-1; h++)
         for(j=1; j<=nlstate*2;j++)
           for(theta=1; theta <=npar; theta++)
           trgradg[h][j][theta]=gradg[h][theta][j];
   
   
        for(i=1;i<=nlstate*2;i++)
         for(j=1;j<=nlstate*2;j++)
           varhe[i][j][(int)age] =0.;
   
       for(h=0;h<=nhstepm-1;h++){
         for(k=0;k<=nhstepm-1;k++){
           matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);
           matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]);
           for(i=1;i<=nlstate*2;i++)
             for(j=1;j<=nlstate*2;j++)
               varhe[i][j][(int)age] += doldm[i][j]*hf*hf;
         }
       }
   
        
       /* Computing expectancies */
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++)        for(j=1; j<=nlstate;j++)
         for (h=0, eij[i][j][(int)age]=0; h<=nhstepm; h++){          for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
           eij[i][j][(int)age] +=p3mat[i][j][h];            eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf;
            
   /* 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]);*/
   
         }          }
      
     hf=1;      fprintf(ficreseij,"%3.0f",age );
     if (stepm >= YEARM) hf=stepm/YEARM;      cptj=0;
     fprintf(ficreseij,"%.0f",age );  
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++){        for(j=1; j<=nlstate;j++){
         fprintf(ficreseij," %.4f", hf*eij[i][j][(int)age]);          cptj++;
           fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) );
       }        }
     fprintf(ficreseij,"\n");      fprintf(ficreseij,"\n");
      
       free_matrix(gm,0,nhstepm,1,nlstate*2);
       free_matrix(gp,0,nhstepm,1,nlstate*2);
       free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*2);
       free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar);
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   }    }
     free_vector(xp,1,npar);
     free_matrix(dnewm,1,nlstate*2,1,npar);
     free_matrix(doldm,1,nlstate*2,1,nlstate*2);
     free_ma3x(varhe,1,nlstate*2,1,nlstate*2,(int) bage, (int)fage);
 }  }
   
 /************ Variance ******************/  /************ Variance ******************/
 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)  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)
 {  {
   /* Variance of health expectancies */    /* Variance of health expectancies */
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/    /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
   double **newm;    double **newm;
   double **dnewm,**doldm;    double **dnewm,**doldm;
   int i, j, nhstepm, hstepm, h;    int i, j, nhstepm, hstepm, h, nstepm ;
   int k, cptcode;    int k, cptcode;
   double *xp;    double *xp;
   double **gp, **gm;    double **gp, **gm;
   double ***gradg, ***trgradg;    double ***gradg, ***trgradg;
   double ***p3mat;    double ***p3mat;
   double age,agelim;    double age,agelim, hf;
   int theta;    int theta;
   
    fprintf(ficresvij,"# Covariances of life expectancies\n");     fprintf(ficresvij,"# Covariances of life expectancies\n");
Line 1557  void varevsij(char fileres[], double *** Line 1733  void varevsij(char fileres[], double ***
   dnewm=matrix(1,nlstate,1,npar);    dnewm=matrix(1,nlstate,1,npar);
   doldm=matrix(1,nlstate,1,nlstate);    doldm=matrix(1,nlstate,1,nlstate);
     
   hstepm=1*YEARM; /* Every year of age */    if(estepm < stepm){
   hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */      printf ("Problem %d lower than %d\n",estepm, stepm);
     }
     else  hstepm=estepm;  
     /* 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.
        nhstepm is the number of hstepm from age to agelim
        nstepm is the number of stepm from age to agelin.
        Look at hpijx to understand the reason of that which relies in memory size
        and note for a fixed period like k years */
     /* 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
        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
        results. So we changed our mind and took the option of the best precision.
     */
     hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */
   agelim = AGESUP;    agelim = AGESUP;
   for (age=bage; age<=fage; age ++){ /* If stepm=6 months */    for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
     nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */      nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
     if (stepm >= YEARM) hstepm=1;      nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
     nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */  
     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
     gradg=ma3x(0,nhstepm,1,npar,1,nlstate);      gradg=ma3x(0,nhstepm,1,npar,1,nlstate);
     gp=matrix(0,nhstepm,1,nlstate);      gp=matrix(0,nhstepm,1,nlstate);
Line 1580  void varevsij(char fileres[], double *** Line 1770  void varevsij(char fileres[], double ***
         for(i=1; i<=nlstate;i++)          for(i=1; i<=nlstate;i++)
           prlim[i][i]=probs[(int)age][i][ij];            prlim[i][i]=probs[(int)age][i][ij];
       }        }
         
       for(j=1; j<= nlstate; j++){        for(j=1; j<= nlstate; j++){
         for(h=0; h<=nhstepm; h++){          for(h=0; h<=nhstepm; h++){
           for(i=1, gp[h][j]=0.;i<=nlstate;i++)            for(i=1, gp[h][j]=0.;i<=nlstate;i++)
Line 1592  void varevsij(char fileres[], double *** Line 1782  void varevsij(char fileres[], double ***
         xp[i] = x[i] - (i==theta ?delti[theta]:0);          xp[i] = x[i] - (i==theta ?delti[theta]:0);
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);          hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
    
       if (popbased==1) {        if (popbased==1) {
         for(i=1; i<=nlstate;i++)          for(i=1; i<=nlstate;i++)
           prlim[i][i]=probs[(int)age][i][ij];            prlim[i][i]=probs[(int)age][i][ij];
Line 1618  void varevsij(char fileres[], double *** Line 1808  void varevsij(char fileres[], double ***
         for(theta=1; theta <=npar; theta++)          for(theta=1; theta <=npar; theta++)
           trgradg[h][j][theta]=gradg[h][theta][j];            trgradg[h][j][theta]=gradg[h][theta][j];
   
       hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
     for(i=1;i<=nlstate;i++)      for(i=1;i<=nlstate;i++)
       for(j=1;j<=nlstate;j++)        for(j=1;j<=nlstate;j++)
         vareij[i][j][(int)age] =0.;          vareij[i][j][(int)age] =0.;
   
     for(h=0;h<=nhstepm;h++){      for(h=0;h<=nhstepm;h++){
       for(k=0;k<=nhstepm;k++){        for(k=0;k<=nhstepm;k++){
         matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);          matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);
         matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);          matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);
         for(i=1;i<=nlstate;i++)          for(i=1;i<=nlstate;i++)
           for(j=1;j<=nlstate;j++)            for(j=1;j<=nlstate;j++)
             vareij[i][j][(int)age] += doldm[i][j];              vareij[i][j][(int)age] += doldm[i][j]*hf*hf;
       }        }
     }      }
     h=1;  
     if (stepm >= YEARM) h=stepm/YEARM;  
     fprintf(ficresvij,"%.0f ",age );      fprintf(ficresvij,"%.0f ",age );
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++){        for(j=1; j<=nlstate;j++){
         fprintf(ficresvij," %.4f", h*vareij[i][j][(int)age]);          fprintf(ficresvij," %.4f", vareij[i][j][(int)age]);
       }        }
     fprintf(ficresvij,"\n");      fprintf(ficresvij,"\n");
     free_matrix(gp,0,nhstepm,1,nlstate);      free_matrix(gp,0,nhstepm,1,nlstate);
Line 1644  void varevsij(char fileres[], double *** Line 1835  void varevsij(char fileres[], double ***
     free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);      free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   } /* End age */    } /* End age */
     
   free_vector(xp,1,npar);    free_vector(xp,1,npar);
   free_matrix(doldm,1,nlstate,1,npar);    free_matrix(doldm,1,nlstate,1,npar);
   free_matrix(dnewm,1,nlstate,1,nlstate);    free_matrix(dnewm,1,nlstate,1,nlstate);
Line 1735  void varprevlim(char fileres[], double * Line 1926  void varprevlim(char fileres[], double *
 }  }
   
 /************ Variance of one-step probabilities  ******************/  /************ Variance of one-step probabilities  ******************/
 void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij)  void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)
 {  {
   int i, j;    int i, j, i1, k1, j1, z1;
   int k=0, cptcode;    int k=0, cptcode;
   double **dnewm,**doldm;    double **dnewm,**doldm;
   double *xp;    double *xp;
Line 1752  void varprob(char fileres[], double **ma Line 1943  void varprob(char fileres[], double **ma
   if((ficresprob=fopen(fileresprob,"w"))==NULL) {    if((ficresprob=fopen(fileresprob,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprob);      printf("Problem with resultfile: %s\n", fileresprob);
   }    }
   printf("Computing variance of one-step probabilities: result on file '%s' \n",fileresprob);    printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
     
   fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n");
     fprintf(ficresprob,"# Age");
     for(i=1; i<=nlstate;i++)
       for(j=1; j<=(nlstate+ndeath);j++)
         fprintf(ficresprob," p%1d-%1d (SE)",i,j);
   
   
     fprintf(ficresprob,"\n");
   
   
   xp=vector(1,npar);    xp=vector(1,npar);
   dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);    dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
   doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));    doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));
     
   cov[1]=1;    cov[1]=1;
   for (age=bage; age<=fage; age ++){    j=cptcoveff;
     cov[2]=age;    if (cptcovn<1) {j=1;ncodemax[1]=1;}
     gradg=matrix(1,npar,1,9);    j1=0;
     trgradg=matrix(1,9,1,npar);    for(k1=1; k1<=1;k1++){
     gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));      for(i1=1; i1<=ncodemax[k1];i1++){
     gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));      j1++;
   
       if  (cptcovn>0) {
         fprintf(ficresprob, "\n#********** Variable ");
         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
         fprintf(ficresprob, "**********\n#");
       }
         
     for(theta=1; theta <=npar; theta++){        for (age=bage; age<=fage; age ++){
       for(i=1; i<=npar; i++)          cov[2]=age;
         xp[i] = x[i] + (i==theta ?delti[theta]:0);          for (k=1; k<=cptcovn;k++) {
                  cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];
       pmij(pmmij,cov,ncovmodel,xp,nlstate);           
      
       k=0;  
       for(i=1; i<= (nlstate+ndeath); i++){  
         for(j=1; j<=(nlstate+ndeath);j++){  
            k=k+1;  
           gp[k]=pmmij[i][j];  
         }          }
       }          for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
           for (k=1; k<=cptcovprod;k++)
       for(i=1; i<=npar; i++)            cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
         xp[i] = x[i] - (i==theta ?delti[theta]:0);         
           gradg=matrix(1,npar,1,9);
           trgradg=matrix(1,9,1,npar);
           gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
           gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
         
           for(theta=1; theta <=npar; theta++){
       pmij(pmmij,cov,ncovmodel,xp,nlstate);            for(i=1; i<=npar; i++)
       k=0;              xp[i] = x[i] + (i==theta ?delti[theta]:0);
       for(i=1; i<=(nlstate+ndeath); i++){           
         for(j=1; j<=(nlstate+ndeath);j++){            pmij(pmmij,cov,ncovmodel,xp,nlstate);
           k=k+1;           
           gm[k]=pmmij[i][j];            k=0;
         }            for(i=1; i<= (nlstate+ndeath); i++){
       }              for(j=1; j<=(nlstate+ndeath);j++){
                 k=k+1;
                 gp[k]=pmmij[i][j];
               }
             }
            
             for(i=1; i<=npar; i++)
               xp[i] = x[i] - (i==theta ?delti[theta]:0);
      
             pmij(pmmij,cov,ncovmodel,xp,nlstate);
             k=0;
             for(i=1; i<=(nlstate+ndeath); i++){
               for(j=1; j<=(nlstate+ndeath);j++){
                 k=k+1;
                 gm[k]=pmmij[i][j];
               }
             }
             
        for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++)            for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++)
            gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];                gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];  
     }          }
   
      for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)  
       for(theta=1; theta <=npar; theta++)  
       trgradg[j][theta]=gradg[theta][j];  
    
      matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);  
      matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);  
   
      pmij(pmmij,cov,ncovmodel,x,nlstate);  
   
      k=0;          for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)
      for(i=1; i<=(nlstate+ndeath); i++){            for(theta=1; theta <=npar; theta++)
        for(j=1; j<=(nlstate+ndeath);j++){              trgradg[j][theta]=gradg[theta][j];
          k=k+1;         
          gm[k]=pmmij[i][j];          matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);
           matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);
          
           pmij(pmmij,cov,ncovmodel,x,nlstate);
          
           k=0;
           for(i=1; i<=(nlstate+ndeath); i++){
             for(j=1; j<=(nlstate+ndeath);j++){
               k=k+1;
               gm[k]=pmmij[i][j];
             }
         }          }
      }  
             
      /*printf("\n%d ",(int)age);       /*printf("\n%d ",(int)age);
      for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){       for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
          
   
        printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));         printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
      }*/       }*/
   
   fprintf(ficresprob,"\n%d ",(int)age);          fprintf(ficresprob,"\n%d ",(int)age);
   
   for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){  
     if (i== 2) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);  
 if (i== 4) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);  
   }  
   
           for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)
             fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i]));
    
         }
       }
     free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));      free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
     free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));      free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
     free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);      free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
     free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);      free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
     }
     free_vector(xp,1,npar);
     fclose(ficresprob);
    
 }  }
  free_vector(xp,1,npar);  
 fclose(ficresprob);  
  exit(0);  
 }  
   
 /***********************************************/  
 /**************** Main Program *****************/  
 /***********************************************/  
   
 /*int main(int argc, char *argv[])*/  /******************* Printing html file ***********/
 int main()  void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \
 {   int lastpass, int stepm, int weightopt, char model[],\
    int imx,int jmin, int jmax, double jmeanint,char optionfile[], \
    char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\
    char version[], int popforecast, int estepm ){
     int jj1, k1, i1, cpt;
     FILE *fichtm;
     /*char optionfilehtm[FILENAMELENGTH];*/
   
   int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;    strcpy(optionfilehtm,optionfile);
   double agedeb, agefin,hf;    strcat(optionfilehtm,".htm");
   double agemin=1.e20, agemax=-1.e20;    if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {
       printf("Problem with %s \n",optionfilehtm), exit(0);
     }
   
   double fret;   fprintf(fichtm,"<body> <font size=\"2\">%s </font> <hr size=\"2\" color=\"#EC5E5E\"> \n
   double **xi,tmp,delta;  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n
   \n
   Total number of observations=%d <br>\n
   Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
   <hr  size=\"2\" color=\"#EC5E5E\">
    <ul><li>Outputs files<br>\n
    - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n
    - Gnuplot file name: <a href=\"%s\">%s</a><br>\n
    - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n
    - 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
    - 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);
   
    fprintf(fichtm,"\n
    - Parameter file with estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>\n
     - Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n
    - Variances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n
    - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>\n
    - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
   
    if(popforecast==1) fprintf(fichtm,"\n
    - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n
    - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n
           <br>",fileres,fileres,fileres,fileres);
    else
      fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model);
   fprintf(fichtm," <li>Graphs</li><p>");
   
   double dum; /* Dummy variable */   m=cptcoveff;
   double ***p3mat;   if (cptcovn < 1) {m=1;ncodemax[1]=1;}
   int *indx;  
   char line[MAXLINE], linepar[MAXLINE];   jj1=0;
   char title[MAXLINE];   for(k1=1; k1<=m;k1++){
   char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];     for(i1=1; i1<=ncodemax[k1];i1++){
   char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];         jj1++;
   char filerest[FILENAMELENGTH];         if (cptcovn > 0) {
   char fileregp[FILENAMELENGTH];           fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
   char popfile[FILENAMELENGTH];           for (cpt=1; cpt<=cptcoveff;cpt++)
   char path[80],pathc[80],pathcd[80],pathtot[80],model[20];             fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
   int firstobs=1, lastobs=10;           fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
   int sdeb, sfin; /* Status at beginning and end */         }
   int c,  h , cpt,l;         fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>
   int ju,jl, mi;  <img src=\"pe%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);    
   int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;         for(cpt=1; cpt<nlstate;cpt++){
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;           fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.gif<br>
   int mobilav=0, fprev, lprev ,fprevfore=1, lprevfore=1,nforecast,popforecast=0;  <img src=\"p%s%d%d.gif\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
   int hstepm, nhstepm;         }
   int *popage;      for(cpt=1; cpt<=nlstate;cpt++) {
          fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
   interval) in state (%d): v%s%d%d.gif <br>
   <img src=\"v%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  
        }
        for(cpt=1; cpt<=nlstate;cpt++) {
           fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.gif <br>
   <img src=\"exp%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
        }
        fprintf(fichtm,"\n<br>- Total life expectancy by age and
   health expectancies in states (1) and (2): e%s%d.gif<br>
   <img src=\"e%s%d.gif\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
   fprintf(fichtm,"\n</body>");
      }
      }
   fclose(fichtm);
   }
   
   /******************* Gnuplot file **************/
   void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
   
     int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
   
     strcpy(optionfilegnuplot,optionfilefiname);
     strcat(optionfilegnuplot,".gp.txt");
     if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
       printf("Problem with file %s",optionfilegnuplot);
     }
   
   #ifdef windows
       fprintf(ficgp,"cd \"%s\" \n",pathc);
   #endif
   m=pow(2,cptcoveff);
    
    /* 1eme*/
     for (cpt=1; cpt<= nlstate ; cpt ++) {
      for (k1=1; k1<= m ; k1 ++) {
   
        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1);
   
   for (i=1; i<= nlstate ; i ++) {
     if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
     else fprintf(ficgp," \%%*lf (\%%*lf)");
   }
       fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1);
       for (i=1; i<= nlstate ; i ++) {
     if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
     else fprintf(ficgp," \%%*lf (\%%*lf)");
   }
     fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"vpl%s\" every :::%d::%d u 1:($2-2*$3) \"\%%lf",fileres,k1-1,k1-1);
        for (i=1; i<= nlstate ; i ++) {
     if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
     else fprintf(ficgp," \%%*lf (\%%*lf)");
   }  
        fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));
   
   fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
      }
     }
     /*2 eme*/
   
     for (k1=1; k1<= m ; k1 ++) {
       fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",ageminpar,fage);
      
       for (i=1; i<= nlstate+1 ; i ++) {
         k=2*i;
         fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:2 \"\%%lf",fileres,k1-1,k1-1);
         for (j=1; j<= nlstate+1 ; j ++) {
     if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
     else fprintf(ficgp," \%%*lf (\%%*lf)");
   }  
         if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");
         else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);
       fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",fileres,k1-1,k1-1);
         for (j=1; j<= nlstate+1 ; j ++) {
           if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
           else fprintf(ficgp," \%%*lf (\%%*lf)");
   }  
         fprintf(ficgp,"\" t\"\" w l 0,");
        fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",fileres,k1-1,k1-1);
         for (j=1; j<= nlstate+1 ; j ++) {
     if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
     else fprintf(ficgp," \%%*lf (\%%*lf)");
   }  
         if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");
         else fprintf(ficgp,"\" t\"\" w l 0,");
       }
       fprintf(ficgp,"\nset out \"e%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),k1);
     }
    
     /*3eme*/
   
     for (k1=1; k1<= m ; k1 ++) {
       for (cpt=1; cpt<= nlstate ; cpt ++) {
         k=2+nlstate*(2*cpt-2);
         fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt);
         /*fprintf(ficgp,",\"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);
   
         }
         fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
       }
       }
    
     /* CV preval stat */
       for (k1=1; k1<= m ; k1 ++) {
       for (cpt=1; cpt<nlstate ; cpt ++) {
         k=3;
         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",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);  
         fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
       }
     }  
    
     /* proba elementaires */
      for(i=1,jk=1; i <=nlstate; i++){
       for(k=1; k <=(nlstate+ndeath); k++){
         if (k != i) {
           for(j=1; j <=ncovmodel; j++){
          
             fprintf(ficgp,"p%d=%f ",jk,p[jk]);
             jk++;
             fprintf(ficgp,"\n");
           }
         }
       }
       }
   
       for(jk=1; jk <=m; jk++) {
     fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",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;
          }
        }
      }
      fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);
      }
      
     fclose(ficgp);
   }  /* end gnuplot */
   
   
   /*************** Moving average **************/
   void movingaverage(double agedeb, double fage,double 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 bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;    double ftolpl=FTOL;
Line 1888  int main() Line 2642  int main()
   double **varpl; /* Variances of prevalence limits by age */    double **varpl; /* Variances of prevalence limits by age */
   double *epj, vepp;    double *epj, vepp;
   double kk1, kk2;    double kk1, kk2;
   double *popeffectif,*popcount;    double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
    
   
   char version[80]="Imach version 64b, May 2001, INED-EUROREVES ";    char version[80]="Imach version 0.8a, May 2002, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];    char *alph[]={"a","a","b","c","d","e"}, str[4];
   
   
Line 1898  int main() Line 2653  int main()
 #include <sys/time.h>  #include <sys/time.h>
 #include <time.h>  #include <time.h>
   char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];    char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];
    
   /* long total_usecs;    /* long total_usecs;
   struct timeval start_time, end_time;    struct timeval start_time, end_time;
     
   gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */    gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
   
   
   printf("\nIMACH, Version 0.7");  
   printf("\nEnter the parameter file name: ");  
   
 #ifdef windows  
   scanf("%s",pathtot);  
   getcwd(pathcd, size);    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);    /*cygwin_split_path(pathtot,path,optionfile);
     printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/      printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
   /* cutv(path,optionfile,pathtot,'\\');*/    /* cutv(path,optionfile,pathtot,'\\');*/
   
 split(pathtot, path,optionfile);    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);    chdir(path);
   replace(pathc,path);    replace(pathc,path);
 #endif  
 #ifdef unix  
   scanf("%s",optionfile);  
 #endif  
   
 /*-------- arguments in the command line --------*/  /*-------- arguments in the command line --------*/
   
   strcpy(fileres,"r");    strcpy(fileres,"r");
   strcat(fileres, optionfile);    strcat(fileres, optionfilefiname);
     strcat(fileres,".txt");    /* Other files have txt extension */
   
   /*---------arguments file --------*/    /*---------arguments file --------*/
   
Line 1949  split(pathtot, path,optionfile); Line 2706  split(pathtot, path,optionfile);
   }    }
   ungetc(c,ficpar);    ungetc(c,ficpar);
   
   fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);    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 ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt,model);    printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d 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 ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt,model);    fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d 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);  
    
   fscanf(ficpar,"fprevalence=%d lprevalence=%d pop_based=%d\n",&fprev,&lprev,&popbased);  
   fprintf(ficparo,"fprevalence=%d lprevalence=%d pop_based=%d\n",fprev,lprev,popbased);  
    
  while((c=getc(ficpar))=='#' && c!= EOF){  
     ungetc(c,ficpar);  
     fgets(line, MAXLINE, ficpar);  
     puts(line);  
     fputs(line,ficparo);  
   }  
   ungetc(c,ficpar);  
    
   fscanf(ficpar,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",&fprevfore,&lprevfore,&nforecast,&mobilav);  
   fprintf(ficparo,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",fprevfore,lprevfore,nforecast,mobilav);  
        
    
 while((c=getc(ficpar))=='#' && c!= EOF){  while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      fgets(line, MAXLINE, ficpar);
Line 1983  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2717  while((c=getc(ficpar))=='#' && c!= EOF){
   }    }
   ungetc(c,ficpar);    ungetc(c,ficpar);
     
   fscanf(ficpar,"popforecast=%d popfile=%s\n",&popforecast,popfile);     
    
   covar=matrix(0,NCOVMAX,1,n);    covar=matrix(0,NCOVMAX,1,n);
   cptcovn=0;    cptcovn=0;
   if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;    if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
Line 2080  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2813  while((c=getc(ficpar))=='#' && c!= EOF){
   printf("\n");    printf("\n");
   
   
     /*-------- data file ----------*/      /*-------- Rewriting paramater file ----------*/
     if((ficres =fopen(fileres,"w"))==NULL) {       strcpy(rfileres,"r");    /* "Rparameterfile */
       printf("Problem with resultfile: %s\n", fileres);goto end;       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);      fprintf(ficres,"#%s\n",version);
         
       /*-------- data file ----------*/
     if((fic=fopen(datafile,"r"))==NULL)    {      if((fic=fopen(datafile,"r"))==NULL)    {
       printf("Problem with datafile: %s\n", datafile);goto end;        printf("Problem with datafile: %s\n", datafile);goto end;
     }      }
Line 2127  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2865  while((c=getc(ficpar))=='#' && c!= EOF){
         cutv(stra, strb,line,' '); moisnais[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);          cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra);
         for (j=ncov;j>=1;j--){          for (j=ncovcol;j>=1;j--){
           cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);            cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
         }          }
         num[i]=atol(stra);          num[i]=atol(stra);
Line 2146  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2884  while((c=getc(ficpar))=='#' && c!= EOF){
     if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;      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[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
     if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;      if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
     }      }*/
     for (i=1; i<=imx; i++) 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]));*/     /*  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*/    /* Calculation of the number of parameter from char model*/
   Tvar=ivector(1,15);    Tvar=ivector(1,15);
   Tprod=ivector(1,15);    Tprod=ivector(1,15);
Line 2163  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2904  while((c=getc(ficpar))=='#' && c!= EOF){
     cptcovn=j+1;      cptcovn=j+1;
     cptcovprod=j1;      cptcovprod=j1;
         
      
     strcpy(modelsav,model);      strcpy(modelsav,model);
     if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){      if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){
       printf("Error. Non available option model=%s ",model);        printf("Error. Non available option model=%s ",model);
Line 2194  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2934  while((c=getc(ficpar))=='#' && c!= EOF){
         }          }
         else {          else {
           cutv(strb,stre,strc,'V');            cutv(strb,stre,strc,'V');
           Tvar[i]=ncov+k1;            Tvar[i]=ncovcol+k1;
           cutv(strb,strc,strd,'V');            cutv(strb,strc,strd,'V');
           Tprod[k1]=i;            Tprod[k1]=i;
           Tvard[k1][1]=atoi(strc);            Tvard[k1][1]=atoi(strc);
Line 2202  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2942  while((c=getc(ficpar))=='#' && c!= EOF){
           Tvar[cptcovn+k2]=Tvard[k1][1];            Tvar[cptcovn+k2]=Tvard[k1][1];
           Tvar[cptcovn+k2+1]=Tvard[k1][2];            Tvar[cptcovn+k2+1]=Tvard[k1][2];
           for (k=1; k<=lastobs;k++)            for (k=1; k<=lastobs;k++)
             covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];              covar[ncovcol+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
           k1++;            k1++;
           k2=k2+2;            k2=k2+2;
         }          }
Line 2219  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2959  while((c=getc(ficpar))=='#' && c!= EOF){
     }      }
 }  }
     
   /*printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);    /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
   printf("cptcovprod=%d ", cptcovprod);    printf("cptcovprod=%d ", cptcovprod);
   scanf("%d ",i);*/    scanf("%d ",i);*/
     fclose(fic);      fclose(fic);
Line 2231  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2971  while((c=getc(ficpar))=='#' && c!= EOF){
     /*-calculation of age at interview from date of interview and age at death -*/      /*-calculation of age at interview from date of interview and age at death -*/
     agev=matrix(1,maxwav,1,imx);      agev=matrix(1,maxwav,1,imx);
   
    for (i=1; i<=imx; i++)      for (i=1; i<=imx; i++) {
      for(m=2; (m<= maxwav); m++)        for(m=2; (m<= maxwav); m++) {
        if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){         if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
          anint[m][i]=9999;           anint[m][i]=9999;
          s[m][i]=-1;           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++)  {      for (i=1; i<=imx; i++)  {
       agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);        agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
       for(m=1; (m<= maxwav); m++){        for(m=1; (m<= maxwav); m++){
         if(s[m][i] >0){          if(s[m][i] >0){
           if (s[m][i] == nlstate+1) {            if (s[m][i] >= nlstate+1) {
             if(agedc[i]>0)              if(agedc[i]>0)
               if(moisdc[i]!=99 && andc[i]!=9999)                if(moisdc[i]!=99 && andc[i]!=9999)
               agev[m][i]=agedc[i];                  agev[m][i]=agedc[i];
             else {              /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
              else {
               if (andc[i]!=9999){                if (andc[i]!=9999){
               printf("Warning negative age at death: %d line:%d\n",num[i],i);                printf("Warning negative age at death: %d line:%d\n",num[i],i);
               agev[m][i]=-1;                agev[m][i]=-1;
Line 2321  printf("Total number of individuals= %d, Line 3065  printf("Total number of individuals= %d,
        for(j=1; j <= ncodemax[k]; j++){         for(j=1; j <= ncodemax[k]; j++){
          for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){           for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){
            h++;             h++;
            if (h>m) h=1;codtab[h][k]=j;             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     /* Calculates basic frequencies. Computes observed prevalence at single age
        and prints on file fileres'p'. */         and prints on file fileres'p'. */
   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprev, lprev,mint,anint);  
       
   free_matrix(mint,1,maxwav,1,n);     
   free_matrix(anint,1,maxwav,1,n);      pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
    
   pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */  
     oldms= 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 */      newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */      savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
Line 2349  printf("Total number of individuals= %d, Line 3101  printf("Total number of individuals= %d,
     }      }
         
     /*--------- results files --------------*/      /*--------- results files --------------*/
     fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, mle,weightopt,model);      fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
    fprintf(ficres,"fprevalence=%d lprevalence=%d pop_based=%d\n",fprev,lprev,popbased);   
    fprintf(ficres,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",fprevfore,lprevfore,nforecast,mobilav);  
   
    jk=1;     jk=1;
    fprintf(ficres,"# Parameters\n");     fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
    printf("# Parameters\n");     printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
    for(i=1,jk=1; i <=nlstate; i++){     for(i=1,jk=1; i <=nlstate; i++){
      for(k=1; k <=(nlstate+ndeath); k++){       for(k=1; k <=(nlstate+ndeath); k++){
        if (k != i)         if (k != i)
Line 2377  printf("Total number of individuals= %d, Line 3128  printf("Total number of individuals= %d,
     ftolhess=ftol; /* Usually correct */      ftolhess=ftol; /* Usually correct */
     hesscov(matcov, p, npar, delti, ftolhess, func);      hesscov(matcov, p, npar, delti, ftolhess, func);
  }   }
     fprintf(ficres,"# Scales\n");      fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
     printf("# Scales\n");      printf("# Scales (for hessian or gradient estimation)\n");
      for(i=1,jk=1; i <=nlstate; i++){       for(i=1,jk=1; i <=nlstate; i++){
       for(j=1; j <=nlstate+ndeath; j++){        for(j=1; j <=nlstate+ndeath; j++){
         if (j!=i) {          if (j!=i) {
           fprintf(ficres,"%1d%1d",i,j);            fprintf(ficres,"%1d%1d",i,j);
           printf("%1d%1d",i,j);            printf("%1d%1d",i,j);
           for(k=1; k<=ncovmodel;k++){            for(k=1; k<=ncovmodel;k++){
             printf(" %.5e",delti[jk]);              printf(" %.5e",delti[jk]);
             fprintf(ficres," %.5e",delti[jk]);              fprintf(ficres," %.5e",delti[jk]);
             jk++;              jk++;
           }            }
           printf("\n");            printf("\n");
           fprintf(ficres,"\n");            fprintf(ficres,"\n");
         }  
       }  
       }  
      
     k=1;  
     fprintf(ficres,"# Covariance\n");  
     printf("# Covariance\n");  
     for(i=1;i<=npar;i++){  
       /*  if (k>nlstate) k=1;  
       i1=(i-1)/(ncovmodel*nlstate)+1;  
       fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);  
       printf("%s%d%d",alph[k],i1,tab[i]);*/  
       fprintf(ficres,"%3d",i);  
       printf("%3d",i);  
       for(j=1; j<=i;j++){  
         fprintf(ficres," %.5e",matcov[i][j]);  
         printf(" %.5e",matcov[i][j]);  
       }  
       fprintf(ficres,"\n");  
       printf("\n");  
       k++;  
     }  
      
     while((c=getc(ficpar))=='#' && c!= EOF){  
       ungetc(c,ficpar);  
       fgets(line, MAXLINE, ficpar);  
       puts(line);  
       fputs(line,ficparo);  
     }  
     ungetc(c,ficpar);  
    
     fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemax, &bage, &fage);  
      
     if (fage <= 2) {  
       bage = agemin;  
       fage = agemax;  
     }  
   
     fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");  
     fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);  
   
      
 /*------------ gnuplot -------------*/  
 chdir(pathcd);  
   if((ficgp=fopen("graph.plt","w"))==NULL) {  
     printf("Problem with file graph.gp");goto end;  
   }  
 #ifdef windows  
   fprintf(ficgp,"cd \"%s\" \n",pathc);  
 #endif  
 m=pow(2,cptcoveff);  
    
  /* 1eme*/  
   for (cpt=1; cpt<= nlstate ; cpt ++) {  
    for (k1=1; k1<= m ; k1 ++) {  
   
 #ifdef windows  
     fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",agemin,fage,fileres,k1-1,k1-1);  
 #endif  
 #ifdef unix  
 fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",agemin,fage,fileres);  
 #endif  
   
 for (i=1; i<= nlstate ; i ++) {  
   if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");  
   else fprintf(ficgp," \%%*lf (\%%*lf)");  
 }  
     fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1);  
     for (i=1; i<= nlstate ; i ++) {  
   if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");  
   else fprintf(ficgp," \%%*lf (\%%*lf)");  
 }  
   fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"vpl%s\" every :::%d::%d u 1:($2-2*$3) \"\%%lf",fileres,k1-1,k1-1);  
      for (i=1; i<= nlstate ; i ++) {  
   if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");  
   else fprintf(ficgp," \%%*lf (\%%*lf)");  
 }    
      fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1));  
 #ifdef unix  
 fprintf(ficgp,"\nset ter gif small size 400,300");  
 #endif  
 fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);  
    }  
   }  
   /*2 eme*/  
   
   for (k1=1; k1<= m ; k1 ++) {  
     fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",agemin,fage);  
      
     for (i=1; i<= nlstate+1 ; i ++) {  
       k=2*i;  
       fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:2 \"\%%lf",fileres,k1-1,k1-1);  
       for (j=1; j<= nlstate+1 ; j ++) {  
   if (j==i) fprintf(ficgp," \%%lf (\%%lf)");  
   else fprintf(ficgp," \%%*lf (\%%*lf)");  
 }    
       if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");  
       else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);  
     fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",fileres,k1-1,k1-1);  
       for (j=1; j<= nlstate+1 ; j ++) {  
         if (j==i) fprintf(ficgp," \%%lf (\%%lf)");  
         else fprintf(ficgp," \%%*lf (\%%*lf)");  
 }    
       fprintf(ficgp,"\" t\"\" w l 0,");  
      fprintf(ficgp,"\"t%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",fileres,k1-1,k1-1);  
       for (j=1; j<= nlstate+1 ; j ++) {  
   if (j==i) fprintf(ficgp," \%%lf (\%%lf)");  
   else fprintf(ficgp," \%%*lf (\%%*lf)");  
 }    
       if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");  
       else fprintf(ficgp,"\" t\"\" w l 0,");  
     }  
     fprintf(ficgp,"\nset out \"e%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),k1);  
   }  
    
   /*3eme*/  
   
   for (k1=1; k1<= m ; k1 ++) {  
     for (cpt=1; cpt<= nlstate ; cpt ++) {  
       k=2+nlstate*(cpt-1);  
       fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",agemin,fage,fileres,k1-1,k1-1,k,cpt);  
       for (i=1; i< nlstate ; i ++) {  
         fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+i,cpt,i+1);  
       }  
       fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);  
     }  
   }  
    
   /* CV preval stat */  
   for (k1=1; k1<= m ; k1 ++) {  
     for (cpt=1; cpt<nlstate ; cpt ++) {  
       k=3;  
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",agemin,agemax,fileres,k1,k+cpt+1,k+1);  
       for (i=1; i< nlstate ; i ++)  
         fprintf(ficgp,"+$%d",k+i+1);  
       fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);  
        
       l=3+(nlstate+ndeath)*cpt;  
       fprintf(ficgp,",\"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",fileres,k1,l+cpt+1,l+1);  
       for (i=1; i< nlstate ; i ++) {  
         l=3+(nlstate+ndeath)*cpt;  
         fprintf(ficgp,"+$%d",l+i+1);  
       }  
       fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);    
       fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);  
     }  
   }    
   
   /* proba elementaires */  
    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,"%s%1d%1d=%f ",alph[j],i,k,p[jk]);*/  
           /*fprintf(ficgp,"%s",alph[1]);*/  
           fprintf(ficgp,"p%d=%f ",jk,p[jk]);  
           jk++;  
           fprintf(ficgp,"\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);
   for(jk=1; jk <=m; jk++) {      estepm=0;
   fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",agemin,agemax);      fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
    i=1;      if (estepm==0 || estepm < stepm) estepm=stepm;
    for(k2=1; k2<=nlstate; k2++) {      if (fage <= 2) {
      k3=i;        bage = ageminpar;
      for(k=1; k<=(nlstate+ndeath); k++) {        fage = agemaxpar;
        if (k != k2){      }
         fprintf(ficgp," exp(p%d+p%d*x",i,i+1);     
 ij=1;      fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
         for(j=3; j <=ncovmodel; j++) {      fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
           if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {      fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
             fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);   
             ij++;      while((c=getc(ficpar))=='#' && c!= EOF){
           }      ungetc(c,ficpar);
           else      fgets(line, MAXLINE, ficpar);
           fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);      puts(line);
         }      fputs(line,ficparo);
           fprintf(ficgp,")/(1");  
          
         for(k1=1; k1 <=nlstate; k1++){    
           fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);  
 ij=1;  
           for(j=3; j <=ncovmodel; j++){  
           if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {  
             fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);  
             ij++;  
           }  
           else  
             fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);  
           }  
           fprintf(ficgp,")");  
         }  
         fprintf(ficgp,") t \"p%d%d\" ", k2,k);  
         if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");  
         i=i+ncovmodel;  
        }  
      }  
    }  
    fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);  
   }    }
        ungetc(c,ficpar);
   fclose(ficgp);  
      
 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);  
     /*  }*/  
      
    /*________fin mle=1_________*/  
      
   
     
     /* No more information from the sample is required now */    fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2);
   /* Reads comments: lines beginning with '#' */    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){    while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      fgets(line, MAXLINE, ficpar);
Line 2622  chdir(path); Line 3203  chdir(path);
     fputs(line,ficparo);      fputs(line,ficparo);
   }    }
   ungetc(c,ficpar);    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);  
     
   fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemax, &bage, &fage);    while((c=getc(ficpar))=='#' && c!= EOF){
   printf("agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax, bage, fage);      ungetc(c,ficpar);
   fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);      fgets(line, MAXLINE, ficpar);
 /*--------- index.htm --------*/      puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
   
   strcpy(optionfilehtm,optionfile);    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);
   strcat(optionfilehtm,".htm");  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);
   if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {  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);
     printf("Problem with %s \n",optionfilehtm);goto end;  
   
   while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
   }    }
     ungetc(c,ficpar);
   
  fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.7 </font> <hr size=\"2\" color=\"#EC5E5E\">    fscanf(ficpar,"popforecast=%d popfile=%s popfiledate=%lf/%lf/%lf last-popfiledate=%lf/%lf/%lf\n",&popforecast,popfile,&jpyram,&mpyram,&anpyram,&jpyram1,&mpyram1,&anpyram1);
 Titre=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>    fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
 Total number of observations=%d <br>    fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
 Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>  
 <hr  size=\"2\" color=\"#EC5E5E\">  
 <li>Outputs files<br><br>\n  
         - Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n  
 - Estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>  
         - Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>  
         - Transition probabilities: <a href=\"pij%s\">pij%s</a><br>  
         - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>  
         - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>  
         - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>  
         - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>  
         - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>  
         - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>  
 <br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);  
   
  fprintf(fichtm," <li>Graphs</li><p>");   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
   
  m=cptcoveff;  /*------------ gnuplot -------------*/
  if (cptcovn < 1) {m=1;ncodemax[1]=1;}   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);
   
  j1=0;  /*--------- index.htm --------*/
  for(k1=1; k1<=m;k1++){  
    for(i1=1; i1<=ncodemax[k1];i1++){    printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast,estepm);
        j1++;  
        if (cptcovn > 0) {  
          fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");  
          for (cpt=1; cpt<=cptcoveff;cpt++)  
            fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[j1][cpt]]);  
          fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");  
        }  
        fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>  
 <img src=\"pe%s%d.gif\">",strtok(optionfile, "."),j1,strtok(optionfile, "."),j1);      
        for(cpt=1; cpt<nlstate;cpt++){  
          fprintf(fichtm,"<br>- Prevalence of disability : p%s%d%d.gif<br>  
 <img src=\"p%s%d%d.gif\">",strtok(optionfile, "."),cpt,j1,strtok(optionfile, "."),cpt,j1);  
        }  
     for(cpt=1; cpt<=nlstate;cpt++) {  
        fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident  
 interval) in state (%d): v%s%d%d.gif <br>  
 <img src=\"v%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,j1,strtok(optionfile, "."),cpt,j1);    
      }  
      for(cpt=1; cpt<=nlstate;cpt++) {  
         fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.gif <br>  
 <img src=\"exp%s%d%d.gif\">",cpt,strtok(optionfile, "."),cpt,j1,strtok(optionfile, "."),cpt,j1);  
      }  
      fprintf(fichtm,"\n<br>- Total life expectancy by age and  
 health expectancies in states (1) and (2): e%s%d.gif<br>  
 <img src=\"e%s%d.gif\">",strtok(optionfile, "."),j1,strtok(optionfile, "."),j1);  
 fprintf(fichtm,"\n</body>");  
    }  
  }  
 fclose(fichtm);  
   
    
   /*--------------- Prevalence limit --------------*/    /*--------------- Prevalence limit --------------*/
     
   strcpy(filerespl,"pl");    strcpy(filerespl,"pl");
Line 2710  fclose(fichtm); Line 3279  fclose(fichtm);
   savms= 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 */    oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
   k=0;    k=0;
   agebase=agemin;    agebase=ageminpar;
   agelim=agemax;    agelim=agemaxpar;
   ftolpl=1.e-10;    ftolpl=1.e-10;
   i1=cptcoveff;    i1=cptcoveff;
   if (cptcovn < 1){i1=1;}    if (cptcovn < 1){i1=1;}
Line 2771  fclose(fichtm); Line 3340  fclose(fichtm);
             for(j=1; j<=nlstate+ndeath;j++)              for(j=1; j<=nlstate+ndeath;j++)
               fprintf(ficrespij," %1d-%1d",i,j);                fprintf(ficrespij," %1d-%1d",i,j);
           fprintf(ficrespij,"\n");            fprintf(ficrespij,"\n");
           for (h=0; h<=nhstepm; h++){             for (h=0; h<=nhstepm; h++){
             fprintf(ficrespij,"%d %.0f %.0f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );              fprintf(ficrespij,"%d %.0f %.0f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
             for(i=1; i<=nlstate;i++)              for(i=1; i<=nlstate;i++)
               for(j=1; j<=nlstate+ndeath;j++)                for(j=1; j<=nlstate+ndeath;j++)
                 fprintf(ficrespij," %.5f", p3mat[i][j][h]);                  fprintf(ficrespij," %.5f", p3mat[i][j][h]);
             fprintf(ficrespij,"\n");              fprintf(ficrespij,"\n");
           }               }
           free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);            free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           fprintf(ficrespij,"\n");            fprintf(ficrespij,"\n");
         }          }
     }      }
   }    }
   
   /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/    varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax);
   
   fclose(ficrespij);    fclose(ficrespij);
   
   /*---------- Forecasting ------------------*/  
   
   strcpy(fileresf,"f");  
   strcat(fileresf,fileres);  
   if((ficresf=fopen(fileresf,"w"))==NULL) {  
     printf("Problem with forecast resultfile: %s\n", fileresf);goto end;  
   }  
   printf("Computing forecasting: result on file '%s' \n", fileresf);  
   
   prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprevfore, lprevfore);  
   
  free_matrix(agev,1,maxwav,1,imx);    /*---------- Forecasting ------------------*/
   /* Mobile average */    if((stepm == 1) && (strcmp(model,".")==0)){
       prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);
   if (cptcoveff==0) ncodemax[cptcoveff]=1;      if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);
     }
   if (mobilav==1) {    else{
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);      erreur=108;
     for (agedeb=bage+3; agedeb<=fage-2; agedeb++)      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);
       for (i=1; i<=nlstate;i++)  
         for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)  
           mobaverage[(int)agedeb][i][cptcod]=0.;  
      
     for (agedeb=bage+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;  
         }  
       }  
     }    
   }    }
   
   stepsize=(int) (stepm+YEARM-1)/YEARM;  
   if (stepm<=12) stepsize=1;  
   
   agelim=AGESUP;  
   hstepm=stepsize*YEARM; /* Every year of age */  
   hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */  
    
   if (popforecast==1) {  
     if((ficpop=fopen(popfile,"r"))==NULL)    {  
       printf("Problem with population file : %s\n",popfile);goto end;  
     }  
     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<=i1;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 Horizon(in years)");  
       for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);  
       if (popforecast==1)  fprintf(ficresf," [Population]");  
   
       for (agedeb=fage; agedeb>=bage; agedeb--){  
         fprintf(ficresf,"\n%.f %.f 0",agedeb, agedeb);  
        if (mobilav==1) {  
         for(j=1; j<=nlstate;j++)  
           fprintf(ficresf," %.3f",mobaverage[(int)agedeb][j][cptcod]);  
         }  
         else {  
           for(j=1; j<=nlstate;j++)  
           fprintf(ficresf," %.3f",probs[(int)agedeb][j][cptcod]);  
         }    
   
        for(j=1; j<=ndeath;j++) fprintf(ficresf," 0.00000");  
        if (popforecast==1) fprintf(ficresf," [%.f] ",popeffectif[(int)agedeb]);  
       }  
        
       for (cpt=1; cpt<=nforecast;cpt++) {  
         fprintf(ficresf,"\n");  
       for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */  
         nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);  
         nhstepm = nhstepm/hstepm;  
         /*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,stepm,hstepm,nhstepm);*/  
   
         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*hstepm/YEARM*stepm==cpt)  
             fprintf(ficresf,"\n%.f %.f %.f",agedeb, agedeb+ h*hstepm/YEARM*stepm, 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][i][cptcod];  
              else kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];  
              if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb];  
             }  
            if (h*hstepm/YEARM*stepm==cpt) {  
              fprintf(ficresf," %.3f", kk1);  
                if (popforecast==1) fprintf(ficresf," [%.f]", kk2);  
            }  
           }  
         }  
         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_imatrix(s,1,maxwav+1,1,n);  
   free_vector(weight,1,n);  
   fclose(ficresf);  
   /*---------- Health expectancies and variances ------------*/    /*---------- Health expectancies and variances ------------*/
   
   strcpy(filerest,"t");    strcpy(filerest,"t");
Line 2944  fclose(fichtm); Line 3392  fclose(fichtm);
     printf("Problem with variance resultfile: %s\n", fileresv);exit(0);      printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
   }    }
   printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);    printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
     calagedate=-1;
   prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
   
   k=0;    k=0;
   for(cptcov=1;cptcov<=i1;cptcov++){    for(cptcov=1;cptcov<=i1;cptcov++){
Line 2956  fclose(fichtm); Line 3406  fclose(fichtm);
   
       fprintf(ficreseij,"\n#****** ");        fprintf(ficreseij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++)        for(j=1;j<=cptcoveff;j++)
         fprintf(ficreseij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);          fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
       fprintf(ficreseij,"******\n");        fprintf(ficreseij,"******\n");
   
       fprintf(ficresvij,"\n#****** ");        fprintf(ficresvij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++)        for(j=1;j<=cptcoveff;j++)
         fprintf(ficresvij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);          fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
       fprintf(ficresvij,"******\n");        fprintf(ficresvij,"******\n");
   
       eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);        eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
       oldm=oldms;savm=savms;        oldm=oldms;savm=savms;
       evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k);          evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov);  
    
       vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);        vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
       oldm=oldms;savm=savms;        oldm=oldms;savm=savms;
       varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);         varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm);
           
   
    
       fprintf(ficrest,"#Total LEs with variances: e.. (std) ");        fprintf(ficrest,"#Total LEs with variances: e.. (std) ");
       for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);        for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
       fprintf(ficrest,"\n");        fprintf(ficrest,"\n");
          
       hf=1;  
       if (stepm >= YEARM) hf=stepm/YEARM;  
       epj=vector(1,nlstate+1);        epj=vector(1,nlstate+1);
       for(age=bage; age <=fage ;age++){        for(age=bage; age <=fage ;age++){
         prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);          prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
Line 2985  fclose(fichtm); Line 3436  fclose(fichtm);
             prlim[i][i]=probs[(int)age][i][k];              prlim[i][i]=probs[(int)age][i][k];
         }          }
                 
         fprintf(ficrest," %.0f",age);          fprintf(ficrest," %4.0f",age);
         for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){          for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
           for(i=1, epj[j]=0.;i <=nlstate;i++) {            for(i=1, epj[j]=0.;i <=nlstate;i++) {
             epj[j] += prlim[i][i]*hf*eij[i][j][(int)age];              epj[j] += prlim[i][i]*eij[i][j][(int)age];
               /*  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
           }            }
           epj[nlstate+1] +=epj[j];            epj[nlstate+1] +=epj[j];
         }          }
   
         for(i=1, vepp=0.;i <=nlstate;i++)          for(i=1, vepp=0.;i <=nlstate;i++)
           for(j=1;j <=nlstate;j++)            for(j=1;j <=nlstate;j++)
             vepp += vareij[i][j][(int)age];              vepp += vareij[i][j][(int)age];
         fprintf(ficrest," %.2f (%.2f)", epj[nlstate+1],hf*sqrt(vepp));          fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
         for(j=1;j <=nlstate;j++){          for(j=1;j <=nlstate;j++){
           fprintf(ficrest," %.2f (%.2f)", epj[j],hf*sqrt(vareij[j][j][(int)age]));            fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
         }          }
         fprintf(ficrest,"\n");          fprintf(ficrest,"\n");
       }        }
     }      }
   }    }
          free_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(ficreseij);    fclose(ficresvij);
  fclose(ficresvij);  
   fclose(ficrest);    fclose(ficrest);
   fclose(ficpar);    fclose(ficpar);
   free_vector(epj,1,nlstate+1);    free_vector(epj,1,nlstate+1);
   /*  scanf("%d ",i); */   
   
   /*------- Variance limit prevalence------*/      /*------- Variance limit prevalence------*/  
   
 strcpy(fileresvpl,"vpl");    strcpy(fileresvpl,"vpl");
   strcat(fileresvpl,fileres);    strcat(fileresvpl,fileres);
   if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {    if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
     printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);      printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);
Line 3024  strcpy(fileresvpl,"vpl"); Line 3475  strcpy(fileresvpl,"vpl");
   }    }
   printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);    printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);
   
  k=0;    k=0;
  for(cptcov=1;cptcov<=i1;cptcov++){    for(cptcov=1;cptcov<=i1;cptcov++){
    for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){      for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
      k=k+1;        k=k+1;
      fprintf(ficresvpl,"\n#****** ");        fprintf(ficresvpl,"\n#****** ");
      for(j=1;j<=cptcoveff;j++)        for(j=1;j<=cptcoveff;j++)
        fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);          fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
      fprintf(ficresvpl,"******\n");        fprintf(ficresvpl,"******\n");
             
      varpl=matrix(1,nlstate,(int) bage, (int) fage);        varpl=matrix(1,nlstate,(int) bage, (int) fage);
      oldm=oldms;savm=savms;        oldm=oldms;savm=savms;
      varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);       varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);
    }      }
  }   }
   
   fclose(ficresvpl);    fclose(ficresvpl);
Line 3055  strcpy(fileresvpl,"vpl"); Line 3506  strcpy(fileresvpl,"vpl");
     
   free_matrix(matcov,1,npar,1,npar);    free_matrix(matcov,1,npar,1,npar);
   free_vector(delti,1,npar);    free_vector(delti,1,npar);
      free_matrix(agev,1,maxwav,1,imx);
   free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);    free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
   
   printf("End of Imach\n");    if(erreur >0)
       printf("End of Imach with error or warning %d\n",erreur);
     else   printf("End of Imach\n");
   /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */    /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */
     
   /* printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);*/    /* printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);*/
Line 3067  strcpy(fileresvpl,"vpl"); Line 3520  strcpy(fileresvpl,"vpl");
   
   
  end:   end:
 #ifdef windows    /* chdir(pathcd);*/
  chdir(pathcd);   /*system("wgnuplot graph.plt");*/
 #endif   /*system("../gp37mgw/wgnuplot graph.plt");*/
     /*system("cd ../gp37mgw");*/
  system("..\\gp37mgw\\wgnuplot graph.plt");   /* system("..\\gp37mgw\\wgnuplot graph.plt");*/
    strcpy(plotcmd,GNUPLOTPROGRAM);
    strcat(plotcmd," ");
    strcat(plotcmd,optionfilegnuplot);
    system(plotcmd);
   
 #ifdef windows   /*#ifdef windows*/
   while (z[0] != 'q') {    while (z[0] != 'q') {
     chdir(pathcd);      /* chdir(path); */
     printf("\nType e to edit output files, c to start again, and q for exiting: ");      printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");
     scanf("%s",z);      scanf("%s",z);
     if (z[0] == 'c') system("./imach");      if (z[0] == 'c') system("./imach");
     else if (z[0] == 'e') {      else if (z[0] == 'e') system(optionfilehtm);
       chdir(path);      else if (z[0] == 'g') system(plotcmd);
       system(optionfilehtm);  
     }  
     else if (z[0] == 'q') exit(0);      else if (z[0] == 'q') exit(0);
   }    }
 #endif    /*#endif */
 }  }
   
   
Removed from v.1.17  
changed lines
  Added in v.1.41.2.1

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