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

Diff for /imach/src/imach.c between versions 1.136 and 1.141

version 1.136, 2010/04/26 20:30:53 version 1.141, 2014/01/26 02:42:01
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.141  2014/01/26 02:42:01  brouard
     * imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested...
   
     Revision 1.140  2011/09/02 10:37:54  brouard
     Summary: times.h is ok with mingw32 now.
   
     Revision 1.139  2010/06/14 07:50:17  brouard
     After the theft of my laptop, I probably lost some lines of codes which were not uploaded to the CVS tree.
     I remember having already fixed agemin agemax which are pointers now but not cvs saved.
   
     Revision 1.138  2010/04/30 18:19:40  brouard
     *** empty log message ***
   
     Revision 1.137  2010/04/29 18:11:38  brouard
     (Module): Checking covariates for more complex models
     than V1+V2. A lot of change to be done. Unstable.
   
   Revision 1.136  2010/04/26 20:30:53  brouard    Revision 1.136  2010/04/26 20:30:53  brouard
   (Module): merging some libgsl code. Fixing computation    (Module): merging some libgsl code. Fixing computation
   of likelione (using inter/intrapolation if mle = 0) in order to    of likelione (using inter/intrapolation if mle = 0) in order to
Line 377 Line 394
 #include <errno.h>  #include <errno.h>
 extern int errno;  extern int errno;
   
 /* #include <sys/time.h> */  #ifdef LINUX
 #include <time.h>  #include <time.h>
 #include "timeval.h"  #include "timeval.h"
   #else
   #include <sys/time.h>
   #endif
   
 #ifdef GSL  #ifdef GSL
 #include <gsl/gsl_errno.h>  #include <gsl/gsl_errno.h>
Line 389  extern int errno; Line 409  extern int errno;
 /* #include <libintl.h> */  /* #include <libintl.h> */
 /* #define _(String) gettext (String) */  /* #define _(String) gettext (String) */
   
 #define MAXLINE 256  #define MAXLINE 1024 /* Was 256. Overflow with 312 with 2 states and 4 covariates. Should be ok */
   
 #define GNUPLOTPROGRAM "gnuplot"  #define GNUPLOTPROGRAM "gnuplot"
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/  /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
Line 423  extern int errno; Line 443  extern int errno;
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
   
 char version[]="Imach version 0.98m, April 2010, INED-EUROREVES-Institut de longevite ";  char version[]="Imach version 0.98n, January 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Sicentific Research 25293121)";
 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 547  double dateintmean=0; Line 567  double dateintmean=0;
   
 double *weight;  double *weight;
 int **s; /* Status */  int **s; /* Status */
 double *agedc, **covar, idx;  double *agedc;
 int **nbcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff;  double  **covar; /**< covar[i,j], value of jth covariate for individual i,
                     * covar=matrix(0,NCOVMAX,1,n); 
                     * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; */
   double  idx; 
   int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
   int **codtab; /**< codtab=imatrix(1,100,1,10); */
   int **Tvard, *Tprod, cptcovprod, *Tvaraff;
 double *lsurv, *lpop, *tpop;  double *lsurv, *lpop, *tpop;
   
 double ftol=FTOL; /* Tolerance for computing Max Likelihood */  double ftol=FTOL; /* Tolerance for computing Max Likelihood */
Line 622  void replace_back_to_slash(char *s, char Line 648  void replace_back_to_slash(char *s, char
 }  }
   
 char *trimbb(char *out, char *in)  char *trimbb(char *out, char *in)
 { /* Trim multiple blanks in line */  { /* Trim multiple blanks in line but keeps first blanks if line starts with blanks */
   char *s;    char *s;
   s=out;    s=out;
   while (*in != '\0'){    while (*in != '\0'){
     while( *in == ' ' && *(in+1) == ' ' && *(in+1) != '\0'){      while( *in == ' ' && *(in+1) == ' '){ /* && *(in+1) != '\0'){*/
       in++;        in++;
     }      }
     *out++ = *in++;      *out++ = *in++;
Line 635  char *trimbb(char *out, char *in) Line 661  char *trimbb(char *out, char *in)
   return s;    return s;
 }  }
   
   char *cutv(char *blocc, char *alocc, char *in, char occ)
   {
     /* cuts string in into blocc and alocc where blocc ends before last occurence of char 'occ' 
        and alocc starts after last occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2')
        gives blocc="abcdef2ghi" and alocc="j".
        If occ is not found blocc is null and alocc is equal to in. Returns alocc
     */
     char *s, *t;
     t=in;s=in;
     while (*in != '\0'){
       while( *in == occ){
         *blocc++ = *in++;
         s=in;
       }
       *blocc++ = *in++;
     }
     if (s == t) /* occ not found */
       *(blocc-(in-s))='\0';
     else
       *(blocc-(in-s)-1)='\0';
     in=s;
     while ( *in != '\0'){
       *alocc++ = *in++;
     }
   
     *alocc='\0';
     return s;
   }
   
 int nbocc(char *s, char occ)  int nbocc(char *s, char occ)
 {  {
   int i,j=0;    int i,j=0;
Line 647  int nbocc(char *s, char occ) Line 702  int nbocc(char *s, char occ)
   return j;    return j;
 }  }
   
 void cutv(char *u,char *v, char*t, char occ)  /* void cutv(char *u,char *v, char*t, char occ) */
 {  /* { */
   /* cuts string t into u and v where u ends before first occurence of char 'occ'   /*   /\* cuts string t into u and v where u ends before last occurence of char 'occ'  */
      and v starts after first occurence of char 'occ' : ex cutv(u,v,"abcdef2ghi2j",'2')  /*      and v starts after last occurence of char 'occ' : ex cutv(u,v,"abcdef2ghi2j",'2') */
      gives u="abcedf" and v="ghi2j" */  /*      gives u="abcdef2ghi" and v="j" *\/ */
   int i,lg,j,p=0;  /*   int i,lg,j,p=0; */
   i=0;  /*   i=0; */
   for(j=0; j<=strlen(t)-1; j++) {  /*   lg=strlen(t); */
     if((t[j]!= occ) && (t[j+1]== occ)) p=j+1;  /*   for(j=0; j<=lg-1; j++) { */
   }  /*     if((t[j]!= occ) && (t[j+1]== occ)) p=j+1; */
   /*   } */
   
   lg=strlen(t);  /*   for(j=0; j<p; j++) { */
   for(j=0; j<p; j++) {  /*     (u[j] = t[j]); */
     (u[j] = t[j]);  /*   } */
   }  /*      u[p]='\0'; */
      u[p]='\0';  
   
    for(j=0; j<= lg; j++) {  /*    for(j=0; j<= lg; j++) { */
     if (j>=(p+1))(v[j-p-1] = t[j]);  /*     if (j>=(p+1))(v[j-p-1] = t[j]); */
   }  /*   } */
 }  /* } */
   
 /********************** nrerror ********************/  /********************** nrerror ********************/
   
Line 1050  char *asc_diff_time(long time_sec, char Line 1105  char *asc_diff_time(long time_sec, char
   sec_left = (sec_left) %(60*60);    sec_left = (sec_left) %(60*60);
   minutes = (sec_left) /60;    minutes = (sec_left) /60;
   sec_left = (sec_left) % (60);    sec_left = (sec_left) % (60);
   sprintf(ascdiff,"%d day(s) %d hour(s) %d minute(s) %d second(s)",days, hours, minutes, sec_left);      sprintf(ascdiff,"%ld day(s) %ld hour(s) %ld minute(s) %ld second(s)",days, hours, minutes, sec_left);  
   return ascdiff;    return ascdiff;
 }  }
   
Line 1227  double **prevalim(double **prlim, int nl Line 1282  double **prevalim(double **prlim, int nl
   for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){    for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
     newm=savm;      newm=savm;
     /* Covariates have to be included here again */      /* Covariates have to be included here again */
      cov[2]=agefin;      cov[2]=agefin;
         
       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 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]]);*/        /*        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++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];      for (k=1; k<=cptcovage;k++) 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]);*/      /*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;
     oldm=newm;      oldm=newm;
     maxmax=0.;      maxmax=0.;
Line 1268  double **prevalim(double **prlim, int nl Line 1323  double **prevalim(double **prlim, int nl
   
 double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )  double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
 {  {
   double s1, s2;    /* According to parameters values stored in x and the covariate's values stored in cov,
        computes the probability to be observed in state j being in state i by appying the
        model to the ncovmodel covariates (including constant and age).
        lnpijopii=ln(pij/pii)= aij+bij*age+cij*v1+dij*v2+... = sum_nc=1^ncovmodel xij(nc)*cov[nc]
        and, according on how parameters are entered, the position of the coefficient xij(nc) of the
        ncth covariate in the global vector x is given by the formula:
        j<i nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel
        j>=i nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel
        Computes ln(pij/pii) (lnpijopii), deduces pij/pii by exponentiation,
        sums on j different of i to get 1-pii/pii, deduces pii, and then all pij.
        Outputs ps[i][j] the probability to be observed in j being in j according to
        the values of the covariates cov[nc] and corresponding parameter values x[nc+shiftij]
     */
     double s1, lnpijopii;
   /*double t34;*/    /*double t34;*/
   int i,j,j1, nc, ii, jj;    int i,j,j1, nc, ii, jj;
   
     for(i=1; i<= nlstate; i++){      for(i=1; i<= nlstate; i++){
       for(j=1; j<i;j++){        for(j=1; j<i;j++){
         for (nc=1, s2=0.;nc <=ncovmodel; nc++){          for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
           /*s2 += param[i][j][nc]*cov[nc];*/            /*lnpijopii += param[i][j][nc]*cov[nc];*/
           s2 += x[(i-1)*nlstate*ncovmodel+(j-1)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];            lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
 /*       printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2); */  /*       printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
         }          }
         ps[i][j]=s2;          ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
 /*      printf("s1=%.17e, s2=%.17e\n",s1,s2); */  /*      printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
       }        }
       for(j=i+1; j<=nlstate+ndeath;j++){        for(j=i+1; j<=nlstate+ndeath;j++){
         for (nc=1, s2=0.;nc <=ncovmodel; nc++){          for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
           s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];            /*lnpijopii += 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); */            lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
   /*        printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
         }          }
         ps[i][j]=s2;          ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
       }        }
     }      }
     /*ps[3][2]=1;*/  
           
     for(i=1; i<= nlstate; i++){      for(i=1; i<= nlstate; i++){
       s1=0;        s1=0;
       for(j=1; j<i; j++){        for(j=1; j<i; j++){
         s1+=exp(ps[i][j]);          s1+=exp(ps[i][j]); /* In fact sums pij/pii */
         /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */          /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
       }        }
       for(j=i+1; j<=nlstate+ndeath; j++){        for(j=i+1; j<=nlstate+ndeath; j++){
         s1+=exp(ps[i][j]);          s1+=exp(ps[i][j]); /* In fact sums pij/pii */
         /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */          /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
       }        }
         /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
       ps[i][i]=1./(s1+1.);        ps[i][i]=1./(s1+1.);
         /* Computing other pijs */
       for(j=1; j<i; j++)        for(j=1; j<i; j++)
         ps[i][j]= exp(ps[i][j])*ps[i][i];          ps[i][j]= exp(ps[i][j])*ps[i][i];
       for(j=i+1; j<=nlstate+ndeath; j++)        for(j=i+1; j<=nlstate+ndeath; j++)
Line 1433  double func( double *x) Line 1503  double func( double *x)
   
   if(mle==1){    if(mle==1){
     for (i=1,ipmx=0, sw=0.; i<=imx; i++){      for (i=1,ipmx=0, sw=0.; i<=imx; i++){
         /* Computes the values of the ncovmodel covariates of the model
            depending if the covariates are fixed or variying (age dependent) and stores them in cov[]
            Then computes with function pmij which return a matrix p[i][j] giving the elementary probability
            to be observed in j being in i according to the model.
          */
       for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];        for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
         /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] 
            is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2] 
            has been calculated etc */
       for(mi=1; mi<= wav[i]-1; mi++){        for(mi=1; mi<= wav[i]-1; mi++){
         for (ii=1;ii<=nlstate+ndeath;ii++)          for (ii=1;ii<=nlstate+ndeath;ii++)
           for (j=1;j<=nlstate+ndeath;j++){            for (j=1;j<=nlstate+ndeath;j++){
Line 1444  double func( double *x) Line 1522  double func( double *x)
           newm=savm;            newm=savm;
           cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;            cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
           for (kk=1; kk<=cptcovage;kk++) {            for (kk=1; kk<=cptcovage;kk++) {
             cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];              cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; /* Tage[kk] gives the data-covariate associated with age */
           }            }
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,            out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                        1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));                         1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
Line 1733  double funcone( double *x) Line 1811  double funcone( double *x)
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;        ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
       /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */        /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
       if(globpr){        if(globpr){
         fprintf(ficresilk,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\          fprintf(ficresilk,"%9ld %6d %2d %2d %1d %1d %3d %11.6f %8.4f\
  %11.6f %11.6f %11.6f ", \   %11.6f %11.6f %11.6f ", \
                 num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],                  num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],
                 2*weight[i]*lli,out[s1][s2],savm[s1][s2]);                  2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
Line 2475  void  concatwav(int wav[], int **dh, int Line 2553  void  concatwav(int wav[], int **dh, int
   }    }
   jmean=sum/k;    jmean=sum/k;
   printf("Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, num[ijmin], jmax, num[ijmax], jmean);    printf("Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, num[ijmin], jmax, num[ijmax], jmean);
   fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean);    fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %d) Max=%d (%d) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean);
  }   }
   
 /*********** Tricode ****************************/  /*********** Tricode ****************************/
Line 2503  void tricode(int *Tvar, int **nbcode, in Line 2581  void tricode(int *Tvar, int **nbcode, in
          (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.*/
     }      }
       /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */
     for (i=0; i<=modmaxcovj; i++) { /* i=-1 ? 0 and 1*/      for (i=0; i<=modmaxcovj; i++) { /* i=-1 ? 0 and 1*//* For each modality of model-cov j */
       if( Ndum[i] != 0 )        if( Ndum[i] != 0 )
         ncodemax[j]++;           ncodemax[j]++; 
       /* Number of modalities of the j th covariate. In fact        /* Number of modalities of the j th covariate. In fact
Line 2515  void tricode(int *Tvar, int **nbcode, in Line 2593  void tricode(int *Tvar, int **nbcode, in
     /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */      /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
     ij=1;       ij=1; 
     for (i=1; i<=ncodemax[j]; i++) { /* i= 1 to 2 for dichotomous */      for (i=1; i<=ncodemax[j]; i++) { /* i= 1 to 2 for dichotomous */
       for (k=0; k<= maxncov; k++) { /* k=-1 ? NCOVMAX*/        for (k=0; k<= modmaxcovj; k++) { /* k=-1 ? NCOVMAX*//* maxncov or modmaxcovj */
         if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */          if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */
           nbcode[Tvar[j]][ij]=k;  /* stores the modality in an array nbcode.             nbcode[Tvar[j]][ij]=k;  /* stores the modality in an array nbcode. 
                                      k is a modality. If we have model=V1+V1*sex                                        k is a modality. If we have model=V1+V1*sex 
Line 2523  void tricode(int *Tvar, int **nbcode, in Line 2601  void tricode(int *Tvar, int **nbcode, in
           ij++;            ij++;
         }          }
         if (ij > ncodemax[j]) break;           if (ij > ncodemax[j]) break; 
       }          }  /* end of loop on */
     }       } /* end of loop on modality */ 
   }      } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/  
     
  for (k=0; k< maxncov; k++) Ndum[k]=0;    for (k=0; k< maxncov; k++) Ndum[k]=0;
     
  for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and age */    for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and age */
    /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/     /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
    ij=Tvar[i]; /* Tvar might be -1 if status was unknown */     ij=Tvar[i]; /* Tvar might be -1 if status was unknown */
    Ndum[ij]++;     Ndum[ij]++;
  }   }
   
  ij=1;   ij=1;
  for (i=1; i<= maxncov; i++) {   for (i=1; i<= maxncov; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
    if((Ndum[i]!=0) && (i<=ncovcol)){     if((Ndum[i]!=0) && (i<=ncovcol)){
      Tvaraff[ij]=i; /*For printing */       Tvaraff[ij]=i; /*For printing */
      ij++;       ij++;
Line 3869  plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1 Line 3947  plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1
                fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);                 fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
              else               else
                fprintf(ficgp," exp(p%d+p%d*x",i,i+1);                 fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
              ij=1;               ij=1;/* To be checked else nbcode[0][0] wrong */
              for(j=3; j <=ncovmodel; j++) {               for(j=3; j <=ncovmodel; j++) {
                if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {                 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]]]);                   fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
Line 4533  int readdata(char datafile[], int firsto Line 4611  int readdata(char datafile[], int firsto
       
   
     for (j=maxwav;j>=1;j--){      for (j=maxwav;j>=1;j--){
       cutv(stra, strb,line,' ');         cutv(stra, strb, line, ' '); 
       if(strb[0]=='.') { /* Missing status */        if(strb[0]=='.') { /* Missing status */
         lval=-1;          lval=-1;
       }else{        }else{
Line 4541  int readdata(char datafile[], int firsto Line 4619  int readdata(char datafile[], int firsto
         lval=strtol(strb,&endptr,10);           lval=strtol(strb,&endptr,10); 
       /*        if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/        /*        if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
         if( strb[0]=='\0' || (*endptr != '\0')){          if( strb[0]=='\0' || (*endptr != '\0')){
           printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);            printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);
           fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);            fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
           return 1;            return 1;
         }          }
       }        }
Line 4556  int readdata(char datafile[], int firsto Line 4634  int readdata(char datafile[], int firsto
         month=99;          month=99;
         year=9999;          year=9999;
       }else{        }else{
         printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);          printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);
         fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);fflush(ficlog);          fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);fflush(ficlog);
         return 1;          return 1;
       }        }
       anint[j][i]= (double) year;         anint[j][i]= (double) year; 
Line 4572  int readdata(char datafile[], int firsto Line 4650  int readdata(char datafile[], int firsto
       month=99;        month=99;
       year=9999;        year=9999;
     }else{      }else{
       printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);        printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);
         fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);          fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);
         return 1;          return 1;
     }      }
     andc[i]=(double) year;       andc[i]=(double) year; 
Line 4587  int readdata(char datafile[], int firsto Line 4665  int readdata(char datafile[], int firsto
       month=99;        month=99;
       year=9999;        year=9999;
     }else{      }else{
       printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);        printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);
       fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);        fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);
         return 1;          return 1;
     }      }
     if (year==9999) {      if (year==9999) {
       printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given.  Exiting.\n",strb, linei,i,line);        printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given.  Exiting.\n",strb, linei,i,line);
       fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog);        fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog);
         return 1;          return 1;
   
     }      }
Line 4605  int readdata(char datafile[], int firsto Line 4683  int readdata(char datafile[], int firsto
     errno=0;      errno=0;
     dval=strtod(strb,&endptr);       dval=strtod(strb,&endptr); 
     if( strb[0]=='\0' || (*endptr != '\0')){      if( strb[0]=='\0' || (*endptr != '\0')){
       printf("Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight.  Exiting.\n",dval, i,line,linei);        printf("Error reading data around '%f' at line number %d, \"%s\" for individual %d\nShould be a weight.  Exiting.\n",dval, i,line,linei);
       fprintf(ficlog,"Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight.  Exiting.\n",dval, i,line,linei);        fprintf(ficlog,"Error reading data around '%f' at line number %d, \"%s\" for individual %d\nShould be a weight.  Exiting.\n",dval, i,line,linei);
       fflush(ficlog);        fflush(ficlog);
       return 1;        return 1;
     }      }
Line 4621  int readdata(char datafile[], int firsto Line 4699  int readdata(char datafile[], int firsto
         errno=0;          errno=0;
         lval=strtol(strb,&endptr,10);           lval=strtol(strb,&endptr,10); 
         if( strb[0]=='\0' || (*endptr != '\0')){          if( strb[0]=='\0' || (*endptr != '\0')){
           printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);            printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);
           fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);fflush(ficlog);            fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);fflush(ficlog);
           return 1;            return 1;
         }          }
       }        }
       if(lval <-1 || lval >1){        if(lval <-1 || lval >1){
         printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\n \          printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \   Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
  for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \   for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
  For example, for multinomial values like 1, 2 and 3,\n \   For example, for multinomial values like 1, 2 and 3,\n \
Line 4636  int readdata(char datafile[], int firsto Line 4714  int readdata(char datafile[], int firsto
  and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \   and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
  output of IMaCh is often meaningless.\n \   output of IMaCh is often meaningless.\n \
  Exiting.\n",lval,linei, i,line,j);   Exiting.\n",lval,linei, i,line,j);
         fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\n \          fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \   Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
  for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \   for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
  For example, for multinomial values like 1, 2 and 3,\n \   For example, for multinomial values like 1, 2 and 3,\n \
Line 4693  int decodemodel ( char model[], int last Line 4771  int decodemodel ( char model[], int last
     cptcovprod=j1; /*Number of products  V1*V2 +v3*age = 2 */      cptcovprod=j1; /*Number of products  V1*V2 +v3*age = 2 */
           
     strcpy(modelsav,model);       strcpy(modelsav,model); 
     if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){      if (strstr(model,"AGE") !=0){
       printf("Error. Non available option model=%s ",model);        printf("Error. AGE must be in lower case 'age' model=%s ",model);
       fprintf(ficlog,"Error. Non available option model=%s ",model);fflush(ficlog);        fprintf(ficlog,"Error. AGE must be in lower case model=%s ",model);fflush(ficlog);
         return 1;
       }
       if (strstr(model,"v") !=0){
         printf("Error. 'v' must be in upper case 'V' model=%s ",model);
         fprintf(ficlog,"Error. 'v' must be in upper case model=%s ",model);fflush(ficlog);
       return 1;        return 1;
     }      }
           
     /* This loop fills the array Tvar from the string 'model'.*/      /* This loop fills the array Tvar from the string 'model'.*/
     /* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */      /* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */
     /*    modelsav=V3*age+V2+V1+V4 strb=V3*age stra=V2+V1+V4       /*   modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4  */
         i=1 Tvar[1]=3 Tage[1]=1        /*  k=4 (age*V3) Tvar[k=4]= 3 (from V3) Tage[cptcovage=1]=4 */
         i=2 Tvar[2]=2      /*  k=3 V4 Tvar[k=3]= 4 (from V4) */
         i=3 Tvar[3]=1      /*  k=2 V1 Tvar[k=2]= 1 (from V1) */
         i=4 Tvar[4]= 4      /*  k=1 Tvar[1]=2 (from V2) */
         i=5 Tvar[5]      /*  k=5 Tvar[5] */
       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]]]; */
      */      /*  } */
     for(k=1; k<=(j+1);k++){      /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
       cutv(strb,stra,modelsav,'+'); /* keeps in strb after the first '+'       for(k=cptcovn; k>=1;k--){
                                      modelsav=V3*age+V2+V1+V4 strb=V3*age stra=V2+V1+V4         cutv(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' 
                                        modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 
                                     */                                       */ 
       /* if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav);*/ /* and analyzes it */        if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
       /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/        /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
       /*scanf("%d",i);*/        /*scanf("%d",i);*/
       if (strchr(strb,'*')) {  /* Model includes a product V3*age+V2+V1+V4 strb=V3*age */        if (strchr(strb,'*')) {  /* Model includes a product V2+V1+V4+V3*age strb=V3*age */
         cutv(strd,strc,strb,'*'); /* strd*strc  Vm*Vn: strb=V3*age strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */          cutv(strd,strc,strb,'*'); /* strd*strc  Vm*Vn: strb=V3*age strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
         if (strcmp(strc,"age")==0) { /* Vn*age */          if (strcmp(strc,"age")==0) { /* Vn*age */
           cptcovprod--;            cptcovprod--;
           cutv(strb,stre,strd,'V'); /* stre="V3" */            cutv(strb,stre,strd,'V'); /* stre="V3" */
           Tvar[k]=atoi(stre);  /* V1+V3*age+V2 Tvar[2]=3, and Tvar[3]=2 */            Tvar[k]=atoi(stre);  /* V2+V1+V4+V3*age Tvar[4]=2 ; V1+V2*age Tvar[2]=2 */
           cptcovage++; /* Sums the number of covariates which include age as a product */            cptcovage++; /* Sums the number of covariates which include age as a product */
           Tage[cptcovage]=k;  /* Tage[1] =2 */            Tage[cptcovage]=k;  /* Tage[1] = 4 */
           /*printf("stre=%s ", stre);*/            /*printf("stre=%s ", stre);*/
         }          } else if (strcmp(strd,"age")==0) { /* or age*Vn */
         else if (strcmp(strd,"age")==0) { /* or age*Vn */  
           cptcovprod--;            cptcovprod--;
           cutv(strb,stre,strc,'V');            cutv(strb,stre,strc,'V');
           Tvar[k]=atoi(stre);            Tvar[k]=atoi(stre);
           cptcovage++;            cptcovage++;
           Tage[cptcovage]=k;            Tage[cptcovage]=k;
         }          } else {  /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2  strb=V3*V2*/
         else {  /* Age is not in the model V1+V3*V2+V2  strb=V3*V2*/            /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
           cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/            cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
           Tvar[k]=ncovcol+k1;  /* find 'n' in Vn and stores in Tvar.             Tvar[k]=ncovcol+k1;  /* For model-covariate k tells which data-covariate to use but
                                   If already ncovcol=2 and model=V2*V1 Tvar[1]=2+1 and Tvar[2]=2+2 etc */                                    because this model-covariate is a construction we invent a new column
                                     ncovcol + k1
                                     If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2
                                     Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */
           cutv(strb,strc,strd,'V'); /* strd was Vm, strc is m */            cutv(strb,strc,strd,'V'); /* strd was Vm, strc is m */
           Tprod[k1]=k;  /* Tprod[1]  */            Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */
           Tvard[k1][1]=atoi(strc); /* m*/            Tvard[k1][1]=atoi(strc); /* m 1 for V1*/
           Tvard[k1][2]=atoi(stre); /* n */            Tvard[k1][2]=atoi(stre); /* n 4 for V4*/
           Tvar[cptcovn+k2]=Tvard[k1][1];            Tvar[cptcovn+k2]=Tvard[k1][1]; /* Tvar[(cptcovn=4+k2=1)=5]= 1 (V1) */
           Tvar[cptcovn+k2+1]=Tvard[k1][2];             Tvar[cptcovn+k2+1]=Tvard[k1][2];  /* Tvar[(cptcovn=4+(k2=1)+1)=6]= 4 (V4) */
           for (i=1; i<=lastobs;i++) /* Computes the new covariate which is a product of covar[n][i]* covar[m][i]            for (i=1; i<=lastobs;i++){
                                      and is stored at ncovol+k1 */              /* Computes the new covariate which is a product of
                  covar[n][i]* covar[m][i] and stores it at ncovol+k1 */
             covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];              covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
             }
           k1++;            k1++;
           k2=k2+2;            k2=k2+2;
         }          } /* End age is not in the model */
       }        } /* End if model includes a product */
       else { /* no more sum */        else { /* no more sum */
         /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/          /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
        /*  scanf("%d",i);*/         /*  scanf("%d",i);*/
         cutv(strd,strc,strb,'V');          cutv(strd,strc,strb,'V');
         Tvar[i]=atoi(strc);          Tvar[k]=atoi(strc);
       }        }
       strcpy(modelsav,stra);  /* modelsav=V2+V3*age+V1+V4 strb=V3*age+V1+V4 */         strcpy(modelsav,stra);  /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ 
       /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);        /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
         scanf("%d",i);*/          scanf("%d",i);*/
     } /* end of loop + */      } /* end of loop + */
Line 4773  int decodemodel ( char model[], int last Line 4861  int decodemodel ( char model[], int last
   scanf("%d ",i);*/    scanf("%d ",i);*/
   
   
   return (0);    return (0); /* with covar[new additional covariate if product] and Tage if age */ 
   endread:    endread:
     printf("Exiting decodemodel: ");      printf("Exiting decodemodel: ");
     return (1);      return (1);
Line 4833  calandcheckages(int imx, int maxwav, dou Line 4921  calandcheckages(int imx, int maxwav, dou
           }            }
           else if(agev[m][i] >*agemax){            else if(agev[m][i] >*agemax){
             *agemax=agev[m][i];              *agemax=agev[m][i];
             /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/              printf(" Max anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.2f\n",m,i,anint[m][i], i,annais[i], *agemax);
           }            }
           /*agev[m][i]=anint[m][i]-annais[i];*/            /*agev[m][i]=anint[m][i]-annais[i];*/
           /*     agev[m][i] = age[i]+2*m;*/            /*     agev[m][i] = age[i]+2*m;*/
Line 4867  calandcheckages(int imx, int maxwav, dou Line 4955  calandcheckages(int imx, int maxwav, dou
 }*/  }*/
   
   
   printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);    printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, *agemin, *agemax);
   fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);     fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, *agemin, *agemax); 
   
   return (0);    return (0);
   endread:    endread:
Line 5089  int main(int argc, char *argv[]) Line 5177  int main(int argc, char *argv[])
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      fgets(line, MAXLINE, ficpar);
     numlinepar++;      numlinepar++;
     puts(line);      fputs(line,stdout);
     fputs(line,ficparo);      fputs(line,ficparo);
     fputs(line,ficlog);      fputs(line,ficlog);
   }    }
Line 5105  int main(int argc, char *argv[]) Line 5193  int main(int argc, char *argv[])
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      fgets(line, MAXLINE, ficpar);
     numlinepar++;      numlinepar++;
     puts(line);      fputs(line, stdout);
       //puts(line);
     fputs(line,ficparo);      fputs(line,ficparo);
     fputs(line,ficlog);      fputs(line,ficlog);
   }    }
Line 5158  int main(int argc, char *argv[]) Line 5247  int main(int argc, char *argv[])
       ungetc(c,ficpar);        ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);        fgets(line, MAXLINE, ficpar);
       numlinepar++;        numlinepar++;
       puts(line);        fputs(line,stdout);
       fputs(line,ficparo);        fputs(line,ficparo);
       fputs(line,ficlog);        fputs(line,ficlog);
     }      }
Line 5208  run imach with mle=-1 to get a correct t Line 5297  run imach with mle=-1 to get a correct t
       ungetc(c,ficpar);        ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);        fgets(line, MAXLINE, ficpar);
       numlinepar++;        numlinepar++;
       puts(line);        fputs(line,stdout);
       fputs(line,ficparo);        fputs(line,ficparo);
       fputs(line,ficlog);        fputs(line,ficlog);
     }      }
Line 5249  run imach with mle=-1 to get a correct t Line 5338  run imach with mle=-1 to get a correct t
       ungetc(c,ficpar);        ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);        fgets(line, MAXLINE, ficpar);
       numlinepar++;        numlinepar++;
       puts(line);        fputs(line,stdout);
       fputs(line,ficparo);        fputs(line,ficparo);
       fputs(line,ficlog);        fputs(line,ficlog);
     }      }
Line 5317  run imach with mle=-1 to get a correct t Line 5406  run imach with mle=-1 to get a correct t
   anint=matrix(1,maxwav,1,n);    anint=matrix(1,maxwav,1,n);
   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,8);    ncodemax=ivector(1,8); /* hard coded ? */
   
   /* Reads data from file datafile */    /* Reads data from file datafile */
   if (readdata(datafile, firstobs, lastobs, &imx)==1)    if (readdata(datafile, firstobs, lastobs, &imx)==1)
     goto end;      goto end;
   
   /* Calculation of the number of parameters from char model */    /* Calculation of the number of parameters from char model */
   Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. Stores the number n of the covariates in Vm+Vn at 1 and m at 2 */      /*    modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4 
   Tprod=ivector(1,15);           k=4 (age*V3) Tvar[k=4]= 3 (from V3) Tag[cptcovage=1]=4
           k=3 V4 Tvar[k=3]= 4 (from V4)
           k=2 V1 Tvar[k=2]= 1 (from V1)
           k=1 Tvar[1]=2 (from V2)
       */
     Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */
     /*  V2+V1+V4+age*V3 is a model with 4 covariates (3 plus signs). 
         For each model-covariate stores the data-covariate id. Tvar[1]=2, Tvar[2]=1, Tvar[3]=4, 
         Tvar[4=age*V3] is 3 and 'age' is recorded in Tage.
     */
     /* For model-covariate k tells which data-covariate to use but
       because this model-covariate is a construction we invent a new column
       ncovcol + k1
       If already ncovcol=4 and model=V2+V1+V1*V4+age*V3
       Tvar[3=V1*V4]=4+1 etc */
     Tprod=ivector(1,15); /* Gives the position of a product */
     /* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3
        if  V2+V1+V1*V4+age*V3+V3*V2   TProd[k1=2]=5 (V3*V2)
     */
   Tvaraff=ivector(1,15);     Tvaraff=ivector(1,15); 
   Tvard=imatrix(1,15,1,2);    Tvard=imatrix(1,15,1,2); /* n=Tvard[k1][1]  and m=Tvard[k1][2] gives the couple n,m of the k1 th product Vn*Vm
   Tage=ivector(1,15);                                    * For V3*V2 (in V2+V1+V1*V4+age*V3+V3*V2), V3*V2 position is 2nd. 
                               * Tvard[k1=2][1]=3 (V3) Tvard[k1=2][2]=2(V2) */
     Tage=ivector(1,15); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age
                            4 covariates (3 plus signs)
                            Tage[1=V3*age]= 4; Tage[2=age*V4] = 3
                         */  
   
   if(decodemodel(model, lastobs) == 1)    if(decodemodel(model, lastobs) == 1)
     goto end;      goto end;
   
     if((double)(lastobs-imx)/(double)imx > 1.10){
       nbwarn++;
       printf("Warning: The value of parameter lastobs=%d is big compared to the \n  effective number of cases imx=%d, please adjust, \n  otherwise you are allocating more memory than necessary.\n",lastobs, imx); 
       fprintf(ficlog,"Warning: The value of parameter lastobs=%d is big compared to the \n  effective number of cases imx=%d, please adjust, \n  otherwise you are allocating more memory than necessary.\n",lastobs, imx); 
     }
     /*  if(mle==1){*/      /*  if(mle==1){*/
   if (weightopt != 1) { /* Maximisation without weights*/    if (weightopt != 1) { /* Maximisation without weights. We can have weights different from 1 but want no weight*/
     for(i=1;i<=n;i++) weight[i]=1.0;      for(i=1;i<=imx;i++) weight[i]=1.0; /* changed to imx */
   }    }
   
     /*-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 -*/
Line 5378  run imach with mle=-1 to get a correct t Line 5495  run imach with mle=-1 to get a correct t
       for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate */        for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate */
         for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){  /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */           for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){  /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */ 
           h++;            h++;
           if (h>m) {            if (h>m) 
             h=1;              h=1;
             codtab[h][k]=j;            codtab[h][k]=j;
             codtab[h][Tvar[k]]=j;            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("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]]);
         }           } 
       }        }
Line 5410  run imach with mle=-1 to get a correct t Line 5526  run imach with mle=-1 to get a correct t
   else{    else{
     fprintf(ficgp,"\n# %s\n", version);       fprintf(ficgp,"\n# %s\n", version); 
     fprintf(ficgp,"# %s\n", optionfilegnuplot);       fprintf(ficgp,"# %s\n", optionfilegnuplot); 
     fprintf(ficgp,"set missing 'NaNq'\n");      //fprintf(ficgp,"set missing 'NaNq'\n");
       fprintf(ficgp,"set datafile missing 'NaNq'\n");
   }    }
   /*  fclose(ficgp);*/    /*  fclose(ficgp);*/
   /*--------- index.htm --------*/    /*--------- index.htm --------*/
Line 5849  Interval (in months) between two waves: Line 5966  Interval (in months) between two waves:
     while((c=getc(ficpar))=='#' && c!= EOF){      while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);        ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);        fgets(line, MAXLINE, ficpar);
       puts(line);        fputs(line,stdout);
       fputs(line,ficparo);        fputs(line,ficparo);
     }      }
     ungetc(c,ficpar);      ungetc(c,ficpar);
Line 5869  Interval (in months) between two waves: Line 5986  Interval (in months) between two waves:
     while((c=getc(ficpar))=='#' && c!= EOF){      while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);        ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);        fgets(line, MAXLINE, ficpar);
       puts(line);        fputs(line,stdout);
       fputs(line,ficparo);        fputs(line,ficparo);
     }      }
     ungetc(c,ficpar);      ungetc(c,ficpar);
Line 5883  Interval (in months) between two waves: Line 6000  Interval (in months) between two waves:
     while((c=getc(ficpar))=='#' && c!= EOF){      while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);        ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);        fgets(line, MAXLINE, ficpar);
       puts(line);        fputs(line,stdout);
       fputs(line,ficparo);        fputs(line,ficparo);
     }      }
     ungetc(c,ficpar);      ungetc(c,ficpar);
Line 5899  Interval (in months) between two waves: Line 6016  Interval (in months) between two waves:
     while((c=getc(ficpar))=='#' && c!= EOF){      while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);        ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);        fgets(line, MAXLINE, ficpar);
       puts(line);        fputs(line,stdout);
       fputs(line,ficparo);        fputs(line,ficparo);
     }      }
     ungetc(c,ficpar);      ungetc(c,ficpar);
Line 6276  Interval (in months) between two waves: Line 6393  Interval (in months) between two waves:
     /*---------- End : free ----------------*/      /*---------- End : free ----------------*/
     if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);      if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);      free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
   }  /* mle==-3 arrives here for freeing */    }  /* mle==-3 arrives here for freeing */
  endfree:   endfree:
   free_matrix(prlim,1,nlstate,1,nlstate);      free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */
     free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);      free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
     free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);      free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
     free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);      free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
Line 6320  Interval (in months) between two waves: Line 6436  Interval (in months) between two waves:
   fprintf(ficlog,"Local time at start %s\nLocal time at end   %s\n",strstart, strtend);     fprintf(ficlog,"Local time at start %s\nLocal time at end   %s\n",strstart, strtend); 
   printf("Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));    printf("Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));
   
   printf("Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec);    printf("Total time was %ld Sec.\n", end_time.tv_sec -start_time.tv_sec);
   fprintf(ficlog,"Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));    fprintf(ficlog,"Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));
   fprintf(ficlog,"Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec);    fprintf(ficlog,"Total time was %ld Sec.\n", end_time.tv_sec -start_time.tv_sec);
   /*  printf("Total time was %d uSec.\n", total_usecs);*/    /*  printf("Total time was %d uSec.\n", total_usecs);*/
 /*   if(fileappend(fichtm,optionfilehtm)){ */  /*   if(fileappend(fichtm,optionfilehtm)){ */
   fprintf(fichtm,"<br>Local time at start %s<br>Local time at end   %s<br>\n</body></html>",strstart, strtend);    fprintf(fichtm,"<br>Local time at start %s<br>Local time at end   %s<br>\n</body></html>",strstart, strtend);
Removed from v.1.136  
changed lines
  Added in v.1.141

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