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

Diff for /imach/src/imach.c between versions 1.191 and 1.197

version 1.191, 2015/07/14 10:00:33 version 1.197, 2015/09/01 18:24:39
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.197  2015/09/01 18:24:39  brouard
     *** empty log message ***
   
     Revision 1.196  2015/08/18 23:17:52  brouard
     Summary: 0.98q5
   
     Revision 1.195  2015/08/18 16:28:39  brouard
     Summary: Adding a hack for testing purpose
   
     After reading the title, ftol and model lines, if the comment line has
     a q, starting with #q, the answer at the end of the run is quit. It
     permits to run test files in batch with ctest. The former workaround was
     $ echo q | imach foo.imach
   
     Revision 1.194  2015/08/18 13:32:00  brouard
     Summary:  Adding error when the covariance matrix doesn't contain the exact number of lines required by the model line.
   
     Revision 1.193  2015/08/04 07:17:42  brouard
     Summary: 0.98q4
   
     Revision 1.192  2015/07/16 16:49:02  brouard
     Summary: Fixing some outputs
   
   Revision 1.191  2015/07/14 10:00:33  brouard    Revision 1.191  2015/07/14 10:00:33  brouard
   Summary: Some fixes    Summary: Some fixes
   
Line 605 Line 628
 /* #define DEBUG */  /* #define DEBUG */
 /* #define DEBUGBRENT */  /* #define DEBUGBRENT */
 #define POWELL /* Instead of NLOPT */  #define POWELL /* Instead of NLOPT */
   #define POWELLF1F3 /* Skip test */
 /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */  /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
 /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */  /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
   
Line 673  typedef struct { Line 697  typedef struct {
 #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */  #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */
 #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */  #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
 #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */  #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */
 #define codtabm(h,k)  1 & (h-1) >> (k-1) ;  #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
 #define MAXN 20000  #define MAXN 20000
 #define YEARM 12. /**< Number of months per year */  #define YEARM 12. /**< Number of months per year */
 #define AGESUP 130  #define AGESUP 130
 #define AGEBASE 40  #define AGEBASE 40
   #define AGEOVERFLOW 1.e20
 #define AGEGOMP 10 /**< Minimal age for Gompertz adjustment */  #define AGEGOMP 10 /**< Minimal age for Gompertz adjustment */
 #ifdef _WIN32  #ifdef _WIN32
 #define DIRSEPARATOR '\\'  #define DIRSEPARATOR '\\'
Line 691  typedef struct { Line 716  typedef struct {
   
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
   #include "version.h"
 char version[]="Imach version 0.98q2, April 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";  char version[]=__IMACH_VERSION__;
   char copyright[]="September 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
 char fullversion[]="$Revision$ $Date$";   char fullversion[]="$Revision$ $Date$"; 
 char strstart[80];  char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];  char optionfilext[10], optionfilefiname[FILENAMELENGTH];
Line 827  int estepm; Line 853  int estepm;
   
 int m,nb;  int m,nb;
 long *num;  long *num;
 int firstpass=0, lastpass=4,*cod, *ncodemax, *Tage,*cens;  int firstpass=0, lastpass=4,*cod, *cens;
   int *ncodemax;  /* ncodemax[j]= Number of modalities of the j th
                      covariate for which somebody answered excluding 
                      undefined. Usually 2: 0 and 1. */
   int *ncodemaxwundef;  /* ncodemax[j]= Number of modalities of the j th
                                covariate for which somebody answered including 
                                undefined. Usually 3: -1, 0 and 1. */
 double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;  double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
 double **pmmij, ***probs;  double **pmmij, ***probs;
 double *ageexmed,*agecens;  double *ageexmed,*agecens;
Line 841  double  **covar; /**< covar[j,i], value Line 873  double  **covar; /**< covar[j,i], value
                   * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */                    * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */
 double  idx;   double  idx; 
 int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */  int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
   int *Tage;
 int *Ndum; /** Freq of modality (tricode */  int *Ndum; /** Freq of modality (tricode */
 int **codtab; /**< codtab=imatrix(1,100,1,10); */  int **codtab; /**< codtab=imatrix(1,100,1,10); */
 int **Tvard, *Tprod, cptcovprod, *Tvaraff;  int **Tvard, *Tprod, cptcovprod, *Tvaraff;
Line 1559  void linmin(double p[], double xi[], int Line 1592  void linmin(double p[], double xi[], int
     xicom[j]=xi[j];       xicom[j]=xi[j]; 
   }     } 
   
   axs=0.0;    /* axs=0.0; */
   xxss=1; /* 1 and using scale */    /* xxss=1; /\* 1 and using scale *\/ */
   xxs=1;    xxs=1;
   do{    /* do{ */
     ax=0.;      ax=0.;
     xx= xxs;      xx= xxs;
     mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);  /* Outputs: xtx[j]=pcom[j]+(*xx)*xicom[j]; fx=f(xtx[j]) */      mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);  /* Outputs: xtx[j]=pcom[j]+(*xx)*xicom[j]; fx=f(xtx[j]) */
Line 1572  void linmin(double p[], double xi[], int Line 1605  void linmin(double p[], double xi[], int
     /* Given input ax=axs and xx=xxs, xx might be too far from ax to get a finite f(xx) */      /* Given input ax=axs and xx=xxs, xx might be too far from ax to get a finite f(xx) */
     /* Searches on line, outputs (ax, xx, bx) such that fx < min(fa and fb) */      /* Searches on line, outputs (ax, xx, bx) such that fx < min(fa and fb) */
     /* Find a bracket a,x,b in direction n=xi ie xicom, order may change. Scale is [0:xxs*xi[j]] et non plus  [0:xi[j]]*/      /* Find a bracket a,x,b in direction n=xi ie xicom, order may change. Scale is [0:xxs*xi[j]] et non plus  [0:xi[j]]*/
     if (fx != fx){    /*   if (fx != fx){ */
         xxs=xxs/scale; /* Trying a smaller xx, closer to initial ax=0 */    /*    xxs=xxs/scale; /\* Trying a smaller xx, closer to initial ax=0 *\/ */
         printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n",  axs, xxs, fx,fb, fa, xx, ax, bx);    /*    printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n",  axs, xxs, fx,fb, fa, xx, ax, bx); */
     }    /*   } */
   }while(fx != fx);    /* }while(fx != fx); */
   
 #ifdef DEBUGLINMIN  #ifdef DEBUGLINMIN
   printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n",  ax,xx,bx,fa,fx,fb);    printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n",  ax,xx,bx,fa,fx,fb);
Line 1653  void powell(double p[], double **xi, int Line 1686  void powell(double p[], double **xi, int
     printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);      printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);
     fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);      fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
 /*     fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */  /*     fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */
    for (i=1;i<=n;i++) {      for (i=1;i<=n;i++) {
       printf(" %d %.12f",i, p[i]);        printf(" %d %.12f",i, p[i]);
       fprintf(ficlog," %d %.12lf",i, p[i]);        fprintf(ficlog," %d %.12lf",i, p[i]);
       fprintf(ficrespow," %.12lf", p[i]);        fprintf(ficrespow," %.12lf", p[i]);
Line 1761  void powell(double p[], double **xi, int Line 1794  void powell(double p[], double **xi, int
       free_vector(ptt,1,n);         free_vector(ptt,1,n); 
       free_vector(pt,1,n);         free_vector(pt,1,n); 
       return;         return; 
     }       } /* enough precision */ 
     if (*iter == ITMAX) nrerror("powell exceeding maximum iterations.");       if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); 
     for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */      for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */
       ptt[j]=2.0*p[j]-pt[j];         ptt[j]=2.0*p[j]-pt[j]; 
Line 1769  void powell(double p[], double **xi, int Line 1802  void powell(double p[], double **xi, int
       pt[j]=p[j];         pt[j]=p[j]; 
     }       } 
     fptt=(*func)(ptt); /* f_3 */      fptt=(*func)(ptt); /* f_3 */
   #ifdef POWELLF1F3
   #else
     if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */      if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */
   #endif
       /* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */        /* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */
       /* From x1 (P0) distance of x2 is at h and x3 is 2h */        /* From x1 (P0) distance of x2 is at h and x3 is 2h */
       /* Let f"(x2) be the 2nd derivative equal everywhere.  */        /* Let f"(x2) be the 2nd derivative equal everywhere.  */
Line 1798  void powell(double p[], double **xi, int Line 1834  void powell(double p[], double **xi, int
       if (t < 0.0) { /* Then we use it for new direction */        if (t < 0.0) { /* Then we use it for new direction */
 #else  #else
       if (directest*t < 0.0) { /* Contradiction between both tests */        if (directest*t < 0.0) { /* Contradiction between both tests */
       printf("directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);          printf("directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
       printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);          printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
       fprintf(ficlog,"directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);          fprintf(ficlog,"directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);
       fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);          fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
     }         } 
       if (directest < 0.0) { /* Then we use it for new direction */        if (directest < 0.0) { /* Then we use it for new direction */
 #endif  #endif
 #ifdef DEBUGLINMIN  #ifdef DEBUGLINMIN
Line 1838  void powell(double p[], double **xi, int Line 1874  void powell(double p[], double **xi, int
         printf("\n");          printf("\n");
         fprintf(ficlog,"\n");          fprintf(ficlog,"\n");
 #endif  #endif
       } /* end of t negative */        } /* end of t or directest negative */
   #ifdef POWELLF1F3
   #else
     } /* end if (fptt < fp)  */      } /* end if (fptt < fp)  */
   }   #endif
     } /* loop iteration */ 
 }   } 
   
 /**** Prevalence limit (stable or period prevalence)  ****************/  /**** Prevalence limit (stable or period prevalence)  ****************/
Line 1873  double **prevalim(double **prlim, int nl Line 1912  double **prevalim(double **prlim, int nl
       cov[3]= agefin*agefin;;        cov[3]= agefin*agefin;;
     for (k=1; k<=cptcovn;k++) {      for (k=1; k<=cptcovn;k++) {
       cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];        cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
       /*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/        /*      printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/
     }      }
     /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */      /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
     for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]*cov[2];      for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]*cov[2];
Line 2875  void lubksb(double **a, int n, int *indx Line 2914  void lubksb(double **a, int n, int *indx
   
 void pstamp(FILE *fichier)  void pstamp(FILE *fichier)
 {  {
   fprintf(fichier,"# %s.%s\n#%s\n#%s\n# %s", optionfilefiname,optionfilext,version,fullversion,strstart);    fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart);
 }  }
   
 /************ Frequencies ********************/  /************ Frequencies ********************/
Line 3298  void tricode(int *Tvar, int **nbcode, in Line 3337  void tricode(int *Tvar, int **nbcode, in
   
   cptcoveff=0;     cptcoveff=0; 
     
   for (k=-1; k < maxncov; k++) Ndum[k]=0;  
   for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */    for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */
   
   /* Loop on covariates without age and products */    /* Loop on covariates without age and products */
   for (j=1; j<=(cptcovs); j++) { /* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only */    for (j=1; j<=(cptcovs); j++) { /* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only */
       for (k=-1; k < maxncov; k++) Ndum[k]=0;
     for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the       for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the 
                                modality of this covariate Vj*/                                  modality of this covariate Vj*/ 
       ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i        ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
Line 3325  void tricode(int *Tvar, int **nbcode, in Line 3364  void tricode(int *Tvar, int **nbcode, in
       /* getting the maximum value of the modality of the covariate        /* getting the maximum value of the modality of the covariate
          (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and           (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and
          female is 1, then modmaxcovj=1.*/           female is 1, then modmaxcovj=1.*/
     } /* end for loop on individuals */      } /* end for loop on individuals i */
     printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj);      printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj);
       fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj);
     cptcode=modmaxcovj;      cptcode=modmaxcovj;
     /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */      /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */
    /*for (i=0; i<=cptcode; i++) {*/     /*for (i=0; i<=cptcode; i++) {*/
     for (i=modmincovj;  i<=modmaxcovj; i++) { /* i=-1 ? 0 and 1*//* For each value of the modality of model-cov j */      for (k=modmincovj;  k<=modmaxcovj; k++) { /* k=-1 ? 0 and 1*//* For each value k of the modality of model-cov j */
       printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], i, Ndum[i]);        printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]);
       if( Ndum[i] != 0 ){ /* Counts if nobody answered, empty modality */        fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]);
         ncodemax[j]++;  /* ncodemax[j]= Number of non-null modalities of the j th covariate. */        if( Ndum[k] != 0 ){ /* Counts if nobody answered modality k ie empty modality, we skip it and reorder */
           if( k != -1){
             ncodemax[j]++;  /* ncodemax[j]= Number of modalities of the j th
                                covariate for which somebody answered excluding 
                                undefined. Usually 2: 0 and 1. */
           }
           ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th
                                covariate for which somebody answered including 
                                undefined. Usually 3: -1, 0 and 1. */
       }        }
       /* In fact  ncodemax[j]=2 (dichotom. variables only) but it could be more for        /* In fact  ncodemax[j]=2 (dichotom. variables only) but it could be more for
          historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */           historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
Line 3350  void tricode(int *Tvar, int **nbcode, in Line 3398  void tricode(int *Tvar, int **nbcode, in
        nbcode[Tvar[j]][1]=0;         nbcode[Tvar[j]][1]=0;
        nbcode[Tvar[j]][2]=1;         nbcode[Tvar[j]][2]=1;
        nbcode[Tvar[j]][3]=2;         nbcode[Tvar[j]][3]=2;
          To be continued (not working yet).
     */      */
     ij=1; /* ij is similar to i but can jumps over null modalities */      ij=0; /* ij is similar to i but can jump over null modalities */
     for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 */      for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/
       for (k=0; k<= cptcode; k++) { /* k=-1 ? k=0 to 1 *//* Could be 1 to 4 */          if (Ndum[i] == 0) { /* If nobody responded to this modality k */
         /*recode from 0 */            break;
         if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */          }
           nbcode[Tvar[j]][ij]=k;  /* stores the modality k in an array nbcode.           ij++;
                                      k is a modality. If we have model=V1+V1*sex           nbcode[Tvar[j]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
                                      then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */          cptcode = ij; /* New max modality for covar j */
           ij++;      } /* end of loop on modality i=-1 to 1 or more */
         }        
         if (ij > ncodemax[j]) break;       /*   for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */
       }  /* end of loop on */      /*  /\*recode from 0 *\/ */
     } /* end of loop on modality */       /*                               k is a modality. If we have model=V1+V1*sex  */
       /*                               then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */
       /*                            But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */
       /*  } */
       /*  /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */
       /*  if (ij > ncodemax[j]) { */
       /*    printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]);  */
       /*    fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */
       /*    break; */
       /*  } */
       /*   }  /\* end of loop on modality k *\/ */
   } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/      } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/  
       
  for (k=-1; k< maxncov; k++) Ndum[k]=0;    for (k=-1; k< maxncov; k++) Ndum[k]=0; 
Line 3374  void tricode(int *Tvar, int **nbcode, in Line 3433  void tricode(int *Tvar, int **nbcode, in
    Ndum[ij]++; /* Might be supersed V1 + V1*age */     Ndum[ij]++; /* Might be supersed V1 + V1*age */
  }    } 
   
  ij=1;   ij=0;
  for (i=0; i<=  maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */   for (i=0; i<=  maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
    /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/     /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
    if((Ndum[i]!=0) && (i<=ncovcol)){     if((Ndum[i]!=0) && (i<=ncovcol)){
        ij++;
      /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/       /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
      Tvaraff[ij]=i; /*For printing (unclear) */       Tvaraff[ij]=i; /*For printing (unclear) */
      ij++;     }else{
    }else         /* Tvaraff[ij]=0; */
        Tvaraff[ij]=0;     }
  }   }
  ij--;   /* ij--; */
  cptcoveff=ij; /*Number of total covariates*/   cptcoveff=ij; /*Number of total covariates*/
   
 }  }
Line 4175  void varprob(char optionfilefiname[], do Line 4235  void varprob(char optionfilefiname[], do
   fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");    fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
   fprintf(fichtm,"\n");    fprintf(fichtm,"\n");
   
   fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of pairs of step probabilities (drawings)</a></h4></li>\n",optionfilehtmcov);    fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4></br>this page is important in order to visualize confidence intervals and especially correlation between disability and recovery</li>\n",optionfilehtmcov);
   fprintf(fichtmcov,"\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n\    fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov);
   file %s<br>\n",optionfilehtmcov);    fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \
   fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated\  
 and drawn. It helps understanding how is the covariance between two incidences.\  and drawn. It helps understanding how is the covariance between two incidences.\
  They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");   They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");
   fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \    fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \
Line 4470  fprintf(fichtm," \n<ul><li><b>Graphs</b> Line 4529  fprintf(fichtm," \n<ul><li><b>Graphs</b>
   
  jj1=0;   jj1=0;
  for(k1=1; k1<=m;k1++){   for(k1=1; k1<=m;k1++){
    for(i1=1; i1<=ncodemax[k1];i1++){     /* for(i1=1; i1<=ncodemax[k1];i1++){ */
      jj1++;       jj1++;
      if (cptcovn > 0) {       if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");         fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++)          for (cpt=1; cpt<=cptcoveff;cpt++){ 
          fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);           fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
            printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);fflush(stdout);
          }
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }       }
      /* Pij */       /* Pij */
Line 4494  fprintf(fichtm," \n<ul><li><b>Graphs</b> Line 4555  fprintf(fichtm," \n<ul><li><b>Graphs</b>
         fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \          fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \
 <img src=\"%s%d%d.png\">",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);  <img src=\"%s%d%d.png\">",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
      }       }
    } /* end i1 */     /* } /\* end i1 *\/ */
  }/* End k1 */   }/* End k1 */
  fprintf(fichtm,"</ul>");   fprintf(fichtm,"</ul>");
   
   
  fprintf(fichtm,"\   fprintf(fichtm,"\
 \n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\  \n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\
  - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n", rfileres,rfileres);   - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br> \
    - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file.<br> \
   But because parameters are usually highly correlated (a higher incidence of disability \
   and a higher incidence of recovery can give very close observed transition) it might \
   be very useful to look not only at linear confidence intervals estimated from the \
   variances but at the covariance matrix. And instead of looking at the estimated coefficients \
   (parameters) of the logistic regression, it might be more meaningful to visualize the \
   covariance matrix of the one-step probabilities. \
   See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres);
   
  fprintf(fichtm," - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",   fprintf(fichtm," - Standard deviation of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
          subdirf2(fileres,"prob"),subdirf2(fileres,"prob"));           subdirf2(fileres,"prob"),subdirf2(fileres,"prob"));
  fprintf(fichtm,"\   fprintf(fichtm,"\
  - Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",   - Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
Line 4544  fprintf(fichtm," \n<ul><li><b>Graphs</b> Line 4612  fprintf(fichtm," \n<ul><li><b>Graphs</b>
   
  jj1=0;   jj1=0;
  for(k1=1; k1<=m;k1++){   for(k1=1; k1<=m;k1++){
    for(i1=1; i1<=ncodemax[k1];i1++){     /* for(i1=1; i1<=ncodemax[k1];i1++){ */
      jj1++;       jj1++;
      if (cptcovn > 0) {       if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");         fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
Line 4563  true period expectancies (those weighted Line 4631  true period expectancies (those weighted
  drawn in addition to the population based expectancies computed using\   drawn in addition to the population based expectancies computed using\
  observed and cahotic prevalences: %s%d.png<br>\   observed and cahotic prevalences: %s%d.png<br>\
 <img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);  <img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
    } /* end i1 */     /* } /\* end i1 *\/ */
  }/* End k1 */   }/* End k1 */
  fprintf(fichtm,"</ul>");   fprintf(fichtm,"</ul>");
  fflush(fichtm);   fflush(fichtm);
Line 4757  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 4825  plot [%.f:%.f]  ", ageminpar, agemaxpar)
                  fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);                   fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
              ij=1;/* To be checked else nbcode[0][0] wrong */               ij=1;/* To be checked else nbcode[0][0] wrong */
              for(j=3; j <=ncovmodel-nagesqr; j++) {               for(j=3; j <=ncovmodel-nagesqr; j++) {
                if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /* Bug valgrind */                 /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
                  fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);                 if(ij <=cptcovage) { /* Bug valgrind */
                  ij++;                   if((j-2)==Tage[ij]) { /* Bug valgrind */
                      fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
                      ij++;
                    }
                }                 }
                else                 else
                  fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);                   fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
Line 4774  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 4845  plot [%.f:%.f]  ", ageminpar, agemaxpar)
       
                ij=1;                 ij=1;
                for(j=3; j <=ncovmodel-nagesqr; j++){                 for(j=3; j <=ncovmodel-nagesqr; j++){
                  if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {                   if(ij <=cptcovage) { /* Bug valgrind */
                    fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);                     if((j-2)==Tage[ij]) { /* Bug valgrind */
                    ij++;                       fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
                        ij++;
                      }
                  }                   }
                  else                   else
                    fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][j-2]]);                     fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
Line 5395  int readdata(char datafile[], int firsto Line 5468  int readdata(char datafile[], int firsto
   
   
   if((fic=fopen(datafile,"r"))==NULL)    {    if((fic=fopen(datafile,"r"))==NULL)    {
     printf("Problem while opening datafile: %s\n", datafile);return 1;      printf("Problem while opening datafile: %s\n", datafile);fflush(stdout);
     fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);return 1;      fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);fflush(ficlog);return 1;
   }    }
   
   i=1;    i=1;
Line 5597  int decodemodel ( char model[], int last Line 5670  int decodemodel ( char model[], int last
   if (strlen(model) >1){ /* If there is at least 1 covariate */    if (strlen(model) >1){ /* If there is at least 1 covariate */
     j=0, j1=0, k1=0, k2=-1, ks=0, cptcovn=0;      j=0, j1=0, k1=0, k2=-1, ks=0, cptcovn=0;
     if (strstr(model,"AGE") !=0){      if (strstr(model,"AGE") !=0){
       printf("Error. AGE must be in lower case 'age' model=1+age+%s ",model);        printf("Error. AGE must be in lower case 'age' model=1+age+%s. ",model);
       fprintf(ficlog,"Error. AGE must be in lower case model=1+age+%s ",model);fflush(ficlog);        fprintf(ficlog,"Error. AGE must be in lower case model=1+age+%s. ",model);fflush(ficlog);
       return 1;        return 1;
     }      }
     if (strstr(model,"v") !=0){      if (strstr(model,"v") !=0){
Line 5611  int decodemodel ( char model[], int last Line 5684  int decodemodel ( char model[], int last
       printf(" strpt=%s, model=%s\n",strpt, model);        printf(" strpt=%s, model=%s\n",strpt, model);
       if(strpt != model){        if(strpt != model){
       printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \        printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
  'model=1+age+age*age+V1' or 'model=1+age+age*age+V1+V1*age', please swap as well as \n \   'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
  corresponding column of parameters.\n",model);   corresponding column of parameters.\n",model);
       fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \        fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
  'model=1+age+age*age+V1' or 'model=1+age+age*age+V1+V1*age', please swap as well as \n \   'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
  corresponding column of parameters.\n",model); fflush(ficlog);   corresponding column of parameters.\n",model); fflush(ficlog);
       return 1;        return 1;
     }      }
Line 6222  int main(int argc, char *argv[]) Line 6295  int main(int argc, char *argv[])
   
   int jj, ll, li, lj, lk;    int jj, ll, li, lj, lk;
   int numlinepar=0; /* Current linenumber of parameter file */    int numlinepar=0; /* Current linenumber of parameter file */
     int num_filled;
   int itimes;    int itimes;
   int NDIM=2;    int NDIM=2;
   int vpopbased=0;    int vpopbased=0;
Line 6231  int main(int argc, char *argv[]) Line 6305  int main(int argc, char *argv[])
   /* FILE *ficgp;*/ /*Gnuplot File */    /* FILE *ficgp;*/ /*Gnuplot File */
   struct stat info;    struct stat info;
   double agedeb=0.;    double agedeb=0.;
   double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;  
     double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW;
   
   double fret;    double fret;
   double dum=0.; /* Dummy variable */    double dum=0.; /* Dummy variable */
Line 6239  int main(int argc, char *argv[]) Line 6314  int main(int argc, char *argv[])
   double ***mobaverage;    double ***mobaverage;
   
   char line[MAXLINE];    char line[MAXLINE];
   char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];    char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
   
     char model[MAXLINE], modeltemp[MAXLINE];
   char pathr[MAXLINE], pathimach[MAXLINE];     char pathr[MAXLINE], pathimach[MAXLINE]; 
   char *tok, *val; /* pathtot */    char *tok, *val; /* pathtot */
   int firstobs=1, lastobs=10;    int firstobs=1, lastobs=10;
   int c,  h , cpt;    int c,  h , cpt, c2;
   int jl=0;    int jl=0;
   int i1, j1, jk, stepsize=0;    int i1, j1, jk, stepsize=0;
     int count=0;
   
   int *tab;     int *tab; 
   int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */    int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
   int mobilav=0,popforecast=0;    int mobilav=0,popforecast=0;
Line 6317  int main(int argc, char *argv[]) Line 6396  int main(int argc, char *argv[])
   getcwd(pathcd, size);    getcwd(pathcd, size);
 #endif  #endif
   syscompilerinfo(0);    syscompilerinfo(0);
   printf("\n%s\n%s",version,fullversion);    printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion);
   if(argc <=1){    if(argc <=1){
     printf("\nEnter the parameter file name: ");      printf("\nEnter the parameter file name: ");
     fgets(pathr,FILENAMELENGTH,stdin);      fgets(pathr,FILENAMELENGTH,stdin);
Line 6381  int main(int argc, char *argv[]) Line 6460  int main(int argc, char *argv[])
     goto end;      goto end;
   }    }
   fprintf(ficlog,"Log filename:%s\n",filelog);    fprintf(ficlog,"Log filename:%s\n",filelog);
   fprintf(ficlog,"\n%s\n%s",version,fullversion);    fprintf(ficlog,"Version %s %s",version,fullversion);
   fprintf(ficlog,"\nEnter the parameter file name: \n");    fprintf(ficlog,"\nEnter the parameter file name: \n");
   fprintf(ficlog,"pathimach=%s\npathtot=%s\n\    fprintf(ficlog,"pathimach=%s\npathtot=%s\n\
  path=%s \n\   path=%s \n\
Line 6389  int main(int argc, char *argv[]) Line 6468  int main(int argc, char *argv[])
  optionfilext=%s\n\   optionfilext=%s\n\
  optionfilefiname='%s'\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);   optionfilefiname='%s'\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
   
   syscompilerinfo(0);    syscompilerinfo(1);
   
   printf("Local time (at start):%s",strstart);    printf("Local time (at start):%s",strstart);
   fprintf(ficlog,"Local time (at start): %s",strstart);    fprintf(ficlog,"Local time (at start): %s",strstart);
Line 6425  int main(int argc, char *argv[]) Line 6504  int main(int argc, char *argv[])
   
   /* Reads comments: lines beginning with '#' */    /* Reads comments: lines beginning with '#' */
   numlinepar=0;    numlinepar=0;
   while((c=getc(ficpar))=='#' && c!= EOF){  
     ungetc(c,ficpar);      /* First parameter line */
     fgets(line, MAXLINE, ficpar);    while(fgets(line, MAXLINE, ficpar)) {
       /* If line starts with a # it is a comment */
       if (line[0] == '#') {
         numlinepar++;
         fputs(line,stdout);
         fputs(line,ficparo);
         fputs(line,ficlog);
         continue;
       }else
         break;
     }
     if((num_filled=sscanf(line,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", \
                           title, datafile, &lastobs, &firstpass,&lastpass)) !=EOF){
       if (num_filled != 5) {
         printf("Should be 5 parameters\n");
       }
     numlinepar++;      numlinepar++;
     fputs(line,stdout);      printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass);
     fputs(line,ficparo);    }
     fputs(line,ficlog);    /* Second parameter line */
     while(fgets(line, MAXLINE, ficpar)) {
       /* If line starts with a # it is a comment */
       if (line[0] == '#') {
         numlinepar++;
         fputs(line,stdout);
         fputs(line,ficparo);
         fputs(line,ficlog);
         continue;
       }else
         break;
     }
     if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \
                           &ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){
       if (num_filled != 8) {
         printf("Not 8\n");
       }
       printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt);
   }    }
   ungetc(c,ficpar);  
   
   fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);    /* Third parameter line */
   numlinepar++;    while(fgets(line, MAXLINE, ficpar)) {
   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=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);      /* If line starts with a # it is a comment */
       if (line[0] == '#') {
         numlinepar++;
         fputs(line,stdout);
         fputs(line,ficparo);
         fputs(line,ficlog);
         continue;
       }else
         break;
     }
     if((num_filled=sscanf(line,"model=1+age%[^.\n]\n", model)) !=EOF){
       if (num_filled != 1) {
         printf("ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line);
         fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line);
         model[0]='\0';
         goto end;
       }
       else{
         if (model[0]=='+'){
           for(i=1; i<=strlen(model);i++)
             modeltemp[i-1]=model[i];
         }
         strcpy(model,modeltemp); 
       }
       printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout);
     }
     /* 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 model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */
     /* numlinepar=numlinepar+3; /\* In general *\/ */
     /* 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=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */
   if(model[strlen(model)-1]=='.') /* Suppressing leading dot in the model */    if(model[strlen(model)-1]=='.') /* Suppressing leading dot in the model */
     model[strlen(model)-1]='\0';      model[strlen(model)-1]='\0';
   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=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);    fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
   fprintf(ficlog,"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=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);    fprintf(ficlog,"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=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
   fflush(ficlog);    fflush(ficlog);
   /* if(model[0]=='#'|| model[0]== '\0'){ */    /* if(model[0]=='#'|| model[0]== '\0'){ */
   if(model[0]=='#'){    if(model[0]=='#'){
Line 6457  int main(int argc, char *argv[]) Line 6595  int main(int argc, char *argv[])
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      fgets(line, MAXLINE, ficpar);
     numlinepar++;      numlinepar++;
       if(line[1]=='q'){ /* This #q will quit imach (the answer is q) */
         z[0]=line[1];
       }
       /* printf("****line [1] = %c \n",line[1]); */
     fputs(line, stdout);      fputs(line, stdout);
     //puts(line);      //puts(line);
     fputs(line,ficparo);      fputs(line,ficparo);
Line 6498  int main(int argc, char *argv[]) Line 6640  int main(int argc, char *argv[])
   }    }
   else if(mle==-3) { /* Main Wizard */    else if(mle==-3) { /* Main Wizard */
     prwizard(ncovmodel, nlstate, ndeath, model, ficparo);      prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
     printf(" You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso);      printf(" You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
     fprintf(ficlog," You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso);      fprintf(ficlog," You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
     param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);      param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
     matcov=matrix(1,npar,1,npar);      matcov=matrix(1,npar,1,npar);
   }    }
Line 6523  int main(int argc, char *argv[]) Line 6665  int main(int argc, char *argv[])
         if(jj==i) continue;          if(jj==i) continue;
         j++;          j++;
         fscanf(ficpar,"%1d%1d",&i1,&j1);          fscanf(ficpar,"%1d%1d",&i1,&j1);
         if ((i1 != i) && (j1 != j)){          if ((i1 != i) || (j1 != jj)){
           printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \            printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \
 It might be a problem of design; if ncovcol and the model are correct\n \  It might be a problem of design; if ncovcol and the model are correct\n \
 run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1);  run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1);
Line 6531  run imach with mle=-1 to get a correct t Line 6673  run imach with mle=-1 to get a correct t
         }          }
         fprintf(ficparo,"%1d%1d",i1,j1);          fprintf(ficparo,"%1d%1d",i1,j1);
         if(mle==1)          if(mle==1)
           printf("%1d%1d",i,j);            printf("%1d%1d",i,jj);
         fprintf(ficlog,"%1d%1d",i,j);          fprintf(ficlog,"%1d%1d",i,jj);
         for(k=1; k<=ncovmodel;k++){          for(k=1; k<=ncovmodel;k++){
           fscanf(ficpar," %lf",&param[i][j][k]);            fscanf(ficpar," %lf",&param[i][j][k]);
           if(mle==1){            if(mle==1){
Line 6613  run imach with mle=-1 to get a correct t Line 6755  run imach with mle=-1 to get a correct t
     for(i=1; i <=npar; i++)      for(i=1; i <=npar; i++)
       for(j=1; j <=npar; j++) matcov[i][j]=0.;        for(j=1; j <=npar; j++) matcov[i][j]=0.;
               
       /* Scans npar lines */
     for(i=1; i <=npar; i++){      for(i=1; i <=npar; i++){
       fscanf(ficpar,"%s",str);        count=fscanf(ficpar,"%1d%1d%1d",&i1,&j1,&jk);
         if(count != 3){
           printf("Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\
   This is probably because your covariance matrix doesn't \n  contain exactly %d lines corresponding to your model line '1+age+%s'.\n\
   Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model);
           fprintf(ficlog,"Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\
   This is probably because your covariance matrix doesn't \n  contain exactly %d lines corresponding to your model line '1+age+%s'.\n\
   Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model);
           exit(1);
         }else
       if(mle==1)        if(mle==1)
         printf("%s",str);          printf("%1d%1d%1d",i1,j1,jk);
       fprintf(ficlog,"%s",str);        fprintf(ficlog,"%1d%1d%1d",i1,j1,jk);
       fprintf(ficparo,"%s",str);        fprintf(ficparo,"%1d%1d%1d",i1,j1,jk);
       for(j=1; j <=i; j++){        for(j=1; j <=i; j++){
         fscanf(ficpar," %le",&matcov[i][j]);          fscanf(ficpar," %le",&matcov[i][j]);
         if(mle==1){          if(mle==1){
Line 6634  run imach with mle=-1 to get a correct t Line 6786  run imach with mle=-1 to get a correct t
       fprintf(ficlog,"\n");        fprintf(ficlog,"\n");
       fprintf(ficparo,"\n");        fprintf(ficparo,"\n");
     }      }
       /* End of read covariance matrix npar lines */
     for(i=1; i <=npar; i++)      for(i=1; i <=npar; i++)
       for(j=i+1;j<=npar;j++)        for(j=i+1;j<=npar;j++)
         matcov[i][j]=matcov[j][i];          matcov[i][j]=matcov[j][i];
Line 6673  run imach with mle=-1 to get a correct t Line 6826  run imach with mle=-1 to get a correct t
   s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */     s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */ 
   tab=ivector(1,NCOVMAX);    tab=ivector(1,NCOVMAX);
   ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */    ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */
     ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */
   
   /* Reads data from file datafile */    /* Reads data from file datafile */
   if (readdata(datafile, firstobs, lastobs, &imx)==1)    if (readdata(datafile, firstobs, lastobs, &imx)==1)
Line 6771  run imach with mle=-1 to get a correct t Line 6925  run imach with mle=-1 to get a correct t
     
   m=pow(2,cptcoveff);    m=pow(2,cptcoveff);
     
   for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */  
     for(i=1; i <=pow(2,cptcoveff-k);i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */   
       for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate ncodemax=2*/  
         for(cpt=1; cpt <=pow(2,k-1); cpt++){  /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */   
           h++;  
           if (h>m)   
             h=1;  
           /**< codtab(h,k)  k   = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1            /**< codtab(h,k)  k   = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1
            * For k=4 covariates, h goes from 1 to 2**k             * For k=4 covariates, h goes from 1 to 2**k
            * codtabm(h,k)=  1 & (h-1) >> (k-1) ;             * codtabm(h,k)=  1 & (h-1) >> (k-1) ;
Line 6791  run imach with mle=-1 to get a correct t Line 6938  run imach with mle=-1 to get a correct t
            *     6     2     1     2     1             *     6     2     1     2     1
            *     7 i=4 1     2     2     1             *     7 i=4 1     2     2     1
            *     8     2     2     2     1             *     8     2     2     2     1
            *     9 i=5 1 i=3 1 i=2 1     1             *     9 i=5 1 i=3 1 i=2 1     2
            *    10     2     1     1     1             *    10     2     1     1     2
            *    11 i=6 1     2     1     1             *    11 i=6 1     2     1     2
            *    12     2     2     1     1             *    12     2     2     1     2
            *    13 i=7 1 i=4 1     2     1                 *    13 i=7 1 i=4 1     2     2    
            *    14     2     1     2     1             *    14     2     1     2     2
            *    15 i=8 1     2     2     1             *    15 i=8 1     2     2     2
            *    16     2     2     2     1             *    16     2     2     2     2
            */             */
           codtab[h][k]=j;    for(h=1; h <=100 ;h++){ 
           /* codtab[12][3]=1; */      /* printf("h=%2d ", h); */
           /*codtab[h][Tvar[k]]=j;*/       for(k=1; k <=10; k++){
           printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]);         /* printf("k=%d %d ",k,codtabm(h,k)); */
         }          codtab[h][k]=codtabm(h,k);
       }       }
     }       /* printf("\n"); */
   }     }
     /* for(k=1;k<=cptcoveff; k++){ /\* scans any effective covariate *\/ */
     /*   for(i=1; i <=pow(2,cptcoveff-k);i++){ /\* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 *\/  */
     /*     for(j=1; j <= ncodemax[k]; j++){ /\* For each modality of this covariate ncodemax=2*\/ */
     /*    for(cpt=1; cpt <=pow(2,k-1); cpt++){  /\* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 *\/  */
     /*      h++; */
     /*      if (h>m)  */
     /*        h=1; */
     /*      codtab[h][k]=j; */
     /*      /\* codtab[12][3]=1; *\/ */
     /*      /\*codtab[h][Tvar[k]]=j;*\/ */
     /*      /\* printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]); *\/ */
     /*    }  */
     /*     } */
     /*   } */
     /* }  */
   /* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]);     /* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]); 
      codtab[1][2]=1;codtab[2][2]=2; */       codtab[1][2]=1;codtab[2][2]=2; */
   /* for(i=1; i <=m ;i++){     /* for(i=1; i <=m ;i++){  */
      for(k=1; k <=cptcovn; k++){    /*    for(k=1; k <=cptcovn; k++){ */
        printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff);    /*      printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); */
      }    /*    } */
      printf("\n");    /*    printf("\n"); */
      }    /* } */
      scanf("%d",i);*/    /*   scanf("%d",i);*/
   
  free_ivector(Ndum,-1,NCOVMAX);   free_ivector(Ndum,-1,NCOVMAX);
   
Line 7072  Interval (in months) between two waves: Line 7234  Interval (in months) between two waves:
     }      }
           
     printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp);      printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp);
     for (i=1;i<=NDIM;i++)       for (i=1;i<=NDIM;i++) {
       printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));        printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
         fprintf(ficlog,"%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
       }
     lsurv=vector(1,AGESUP);      lsurv=vector(1,AGESUP);
     lpop=vector(1,AGESUP);      lpop=vector(1,AGESUP);
     tpop=vector(1,AGESUP);      tpop=vector(1,AGESUP);
Line 7106  Interval (in months) between two waves: Line 7269  Interval (in months) between two waves:
           
           
     replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */      replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */
     printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);      if(ageminpar == AGEOVERFLOW ||agemaxpar == AGEOVERFLOW){
               printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\
   This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
   Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
           fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\
   This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
   Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
       }else
         printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
     printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \      printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \
                      stepm, weightopt,\                       stepm, weightopt,\
                      model,imx,p,matcov,agemortsup);                       model,imx,p,matcov,agemortsup);
Line 7142  Interval (in months) between two waves: Line 7312  Interval (in months) between two waves:
     }      }
           
     /*--------- results files --------------*/      /*--------- results files --------------*/
     fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);      fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
           
           
     fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");      fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
Line 7171  Interval (in months) between two waves: Line 7341  Interval (in months) between two waves:
       ftolhess=ftol; /* Usually correct */        ftolhess=ftol; /* Usually correct */
       hesscov(matcov, p, npar, delti, ftolhess, func);        hesscov(matcov, p, npar, delti, ftolhess, func);
     }      }
       printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
       fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n  It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
       for(i=1,jk=1; i <=nlstate; i++){
         for(k=1; k <=(nlstate+ndeath); k++){
           if (k != i) {
             printf("%d%d ",i,k);
             fprintf(ficlog,"%d%d ",i,k);
             for(j=1; j <=ncovmodel; j++){
               printf("%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
               fprintf(ficlog,"%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
               jk++; 
             }
             printf("\n");
             fprintf(ficlog,"\n");
           }
         }
       }
   
     fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");      fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
     printf("# Scales (for hessian or gradient estimation)\n");      printf("# Scales (for hessian or gradient estimation)\n");
     fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");      fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");
Line 7327  Interval (in months) between two waves: Line 7515  Interval (in months) between two waves:
     dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.;      dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.;
           
     fscanf(ficpar,"pop_based=%d\n",&popbased);      fscanf(ficpar,"pop_based=%d\n",&popbased);
       fprintf(ficlog,"pop_based=%d\n",popbased);
     fprintf(ficparo,"pop_based=%d\n",popbased);         fprintf(ficparo,"pop_based=%d\n",popbased);   
     fprintf(ficres,"pop_based=%d\n",popbased);         fprintf(ficres,"pop_based=%d\n",popbased);   
           
Line 7351  Interval (in months) between two waves: Line 7540  Interval (in months) between two waves:
     /* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */      /* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */
           
     replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */      replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */
     printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);      if(ageminpar == AGEOVERFLOW ||agemaxpar == -AGEOVERFLOW){
           printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\
   This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
   Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
           fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\
   This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
   Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
       }else
         printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
           
     printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\      printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\
                  model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\                   model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\
Line 7646  Interval (in months) between two waves: Line 7843  Interval (in months) between two waves:
     free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);      free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
   
     free_ivector(ncodemax,1,NCOVMAX);      free_ivector(ncodemax,1,NCOVMAX);
       free_ivector(ncodemaxwundef,1,NCOVMAX);
     free_ivector(Tvar,1,NCOVMAX);      free_ivector(Tvar,1,NCOVMAX);
     free_ivector(Tprod,1,NCOVMAX);      free_ivector(Tprod,1,NCOVMAX);
     free_ivector(Tvaraff,1,NCOVMAX);      free_ivector(Tvaraff,1,NCOVMAX);
Line 7753  Interval (in months) between two waves: Line 7951  Interval (in months) between two waves:
   }    }
   end:    end:
   while (z[0] != 'q') {    while (z[0] != 'q') {
     printf("\nType  q for exiting: ");      printf("\nType  q for exiting: "); fflush(stdout);
     scanf("%s",z);      scanf("%s",z);
   }    }
 }  }
Removed from v.1.191  
changed lines
  Added in v.1.197

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