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

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

version 1.135, 2009/10/29 15:33:14 version 1.136, 2010/04/26 20:30:53
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.136  2010/04/26 20:30:53  brouard
     (Module): merging some libgsl code. Fixing computation
     of likelione (using inter/intrapolation if mle = 0) in order to
     get same likelihood as if mle=1.
     Some cleaning of code and comments added.
   
   Revision 1.135  2009/10/29 15:33:14  brouard    Revision 1.135  2009/10/29 15:33:14  brouard
   (Module): Now imach stops if date of birth, at least year of birth, is not given. Some cleaning of the code.    (Module): Now imach stops if date of birth, at least year of birth, is not given. Some cleaning of the code.
   
Line 375  extern int errno; Line 381  extern int errno;
 #include <time.h>  #include <time.h>
 #include "timeval.h"  #include "timeval.h"
   
   #ifdef GSL
   #include <gsl/gsl_errno.h>
   #include <gsl/gsl_multimin.h>
   #endif
   
 /* #include <libintl.h> */  /* #include <libintl.h> */
 /* #define _(String) gettext (String) */  /* #define _(String) gettext (String) */
   
Line 412  extern int errno; Line 423  extern int errno;
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
   
 char version[]="Imach version 0.98l, October 2009, INED-EUROREVES-Institut de longevite ";  char version[]="Imach version 0.98m, April 2010, INED-EUROREVES-Institut de longevite ";
 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 439  int **bh; /* bh[mi][i] is the bias (+ or Line 450  int **bh; /* bh[mi][i] is the bias (+ or
 double jmean=1; /* Mean space between 2 waves */  double jmean=1; /* Mean space between 2 waves */
 double **oldm, **newm, **savm; /* Working pointers to matrices */  double **oldm, **newm, **savm; /* Working pointers to matrices */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */  double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
 FILE *fic,*ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;  /*FILE *fic ; */ /* Used in readdata only */
   FILE *ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
 FILE *ficlog, *ficrespow;  FILE *ficlog, *ficrespow;
 int globpr=0; /* Global variable for printing or not */  int globpr=0; /* Global variable for printing or not */
 double fretone; /* Only one call to likelihood */  double fretone; /* Only one call to likelihood */
Line 1223  double **prevalim(double **prlim, int nl Line 1235  double **prevalim(double **prlim, int nl
       }        }
       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]);*/
Line 1712  double funcone( double *x) Line 1724  double funcone( double *x)
         lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */          lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
       } else if (mle==4){  /* mle=4 no inter-extrapolation */        } else if (mle==4){  /* mle=4 no inter-extrapolation */
         lli=log(out[s1][s2]); /* Original formula */          lli=log(out[s1][s2]); /* Original formula */
       } else{  /* ml>=5 no inter-extrapolation no jackson =0.8a */        } else{  /* mle=0 back to 1 */
         lli=log(out[s1][s2]); /* Original formula */          lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
           /*lli=log(out[s1][s2]); */ /* Original formula */
       } /* End of if */        } /* End of if */
       ipmx +=1;        ipmx +=1;
       sw += weight[i];        sw += weight[i];
Line 2436  void  concatwav(int wav[], int **dh, int Line 2449  void  concatwav(int wav[], int **dh, int
             dh[mi][i]=jk;              dh[mi][i]=jk;
             bh[mi][i]=0;              bh[mi][i]=0;
           }else{ /* We want a negative bias in order to only have interpolation ie            }else{ /* We want a negative bias in order to only have interpolation ie
                   * at the price of an extra matrix product in likelihood */                    * to avoid the price of an extra matrix product in likelihood */
             dh[mi][i]=jk+1;              dh[mi][i]=jk+1;
             bh[mi][i]=ju;              bh[mi][i]=ju;
           }            }
Line 2468  void  concatwav(int wav[], int **dh, int Line 2481  void  concatwav(int wav[], int **dh, int
 /*********** Tricode ****************************/  /*********** Tricode ****************************/
 void tricode(int *Tvar, int **nbcode, int imx)  void tricode(int *Tvar, int **nbcode, int imx)
 {  {
       /* Uses cptcovn+2*cptcovprod as the number of covariates */
   /*      Tvar[i]=atoi(stre); /* find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 */    /*      Tvar[i]=atoi(stre); /* find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 */
   
   int Ndum[20],ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;    int Ndum[20],ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
   int cptcode=0;    int modmaxcovj=0; /* Modality max of covariates j */
   cptcoveff=0;     cptcoveff=0; 
     
   for (k=0; k<maxncov; k++) Ndum[k]=0;    for (k=0; k<maxncov; k++) Ndum[k]=0;
   for (k=1; k<=7; k++) ncodemax[k]=0; /* Horrible constant again */    for (k=1; k<=7; k++) ncodemax[k]=0; /* Horrible constant again */
   
   for (j=1; j<=(cptcovn+2*cptcovprod); j++) { /* For each covariate */    for (j=1; j<=(cptcovn+2*cptcovprod); j++) { /* For each covariate j */
     for (i=1; i<=imx; i++) { /*reads the data file to get the maximum       for (i=1; i<=imx; i++) { /*reads the data file to get the maximum value of the 
                                modality*/                                  modality of this covariate Vj*/ 
       ij=(int)(covar[Tvar[j]][i]); /* ij is the modality of this individual, might be -1*/        ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Finds for covariate j, n=Tvar[j] of Vn . ij is the
       Ndum[ij]++; /*counts the occurence of this modality */                                        modality of the nth covariate of individual i. */
         Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/
       /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/        /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
       if (ij > cptcode) cptcode=ij; /* getting the maximum value of the modality of the covariate  (should be 0 or 1 now)         if (ij > modmaxcovj) modmaxcovj=ij; 
                                        Tvar[j]. If V=sex and male is 0 and         /* getting the maximum value of the modality of the covariate
                                        female is 1, then  cptcode=1.*/           (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and
            female is 1, then modmaxcovj=1.*/
     }      }
   
     for (i=0; i<=cptcode; i++) { /* i=-1 ?*/      for (i=0; i<=modmaxcovj; i++) { /* i=-1 ? 0 and 1*/
       if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j        if( Ndum[i] != 0 )
                                        th covariate. In fact          ncodemax[j]++; 
                                        ncodemax[j]=2        /* Number of modalities of the j th covariate. In fact
                                        (dichotom. variables only) but           ncodemax[j]=2 (dichotom. variables only) but it could be more for
                                        it can be more */           historical reasons */
     } /* Ndum[-1] number of undefined modalities */      } /* Ndum[-1] number of undefined modalities */
   
       /* 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 (i=1; i<=ncodemax[j]; i++) { /* i= 1 to 2 for dichotomous */
       for (k=0; k<= maxncov; k++) { /* k=-1 ?*/        for (k=0; k<= maxncov; k++) { /* k=-1 ? NCOVMAX*/
         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 4386  double gompertz(double x[]) Line 4402  double gompertz(double x[])
   return -2*L*num/sump;    return -2*L*num/sump;
 }  }
   
   #ifdef GSL
   /******************* Gompertz_f Likelihood ******************************/
   double gompertz_f(const gsl_vector *v, void *params)
   { 
     double A,B,LL=0.0,sump=0.,num=0.;
     double *x= (double *) v->data;
     int i,n=0; /* n is the size of the sample */
   
     for (i=0;i<=imx-1 ; i++) {
       sump=sump+weight[i];
       /*    sump=sump+1;*/
       num=num+1;
     }
    
    
     /* for (i=0; i<=imx; i++) 
        if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/
     printf("x[0]=%lf x[1]=%lf\n",x[0],x[1]);
     for (i=1;i<=imx ; i++)
       {
         if (cens[i] == 1 && wav[i]>1)
           A=-x[0]/(x[1])*(exp(x[1]*(agecens[i]-agegomp))-exp(x[1]*(ageexmed[i]-agegomp)));
         
         if (cens[i] == 0 && wav[i]>1)
           A=-x[0]/(x[1])*(exp(x[1]*(agedc[i]-agegomp))-exp(x[1]*(ageexmed[i]-agegomp)))
                +log(x[0]/YEARM)+x[1]*(agedc[i]-agegomp)+log(YEARM);  
         
         /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
         if (wav[i] > 1 ) { /* ??? */
           LL=LL+A*weight[i];
           /*      printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/
         }
       }
   
    /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/
     printf("x[0]=%lf x[1]=%lf -2*LL*num/sump=%lf\n",x[0],x[1],-2*LL*num/sump);
    
     return -2*LL*num/sump;
   }
   #endif
   
 /******************* Printing html file ***********/  /******************* Printing html file ***********/
 void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \  void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \
                   int lastpass, int stepm, int weightopt, char model[],\                    int lastpass, int stepm, int weightopt, char model[],\
Line 4433  void printinggnuplotmort(char fileres[], Line 4490  void printinggnuplotmort(char fileres[],
   
 }   } 
   
   int readdata(char datafile[], int firstobs, int lastobs, int *imax)
   {
   
     /*-------- data file ----------*/
     FILE *fic;
     char dummy[]="                         ";
     int i, j, n;
     int linei, month, year,iout;
     char line[MAXLINE], linetmp[MAXLINE];
     char stra[80], strb[80];
     char *stratrunc;
     int lstra;
   
   
     if((fic=fopen(datafile,"r"))==NULL)    {
       printf("Problem while opening datafile: %s\n", datafile);return 1;
       fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);return 1;
     }
   
 /***********************************************/    i=1;
 /**************** Main Program *****************/    linei=0;
 /***********************************************/    while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) {
       linei=linei+1;
       for(j=strlen(line); j>=0;j--){  /* Untabifies line */
         if(line[j] == '\t')
           line[j] = ' ';
       }
       for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10)||(line[j]==13);j--){
         ;
       };
       line[j+1]=0;  /* Trims blanks at end of line */
       if(line[0]=='#'){
         fprintf(ficlog,"Comment line\n%s\n",line);
         printf("Comment line\n%s\n",line);
         continue;
       }
       trimbb(linetmp,line); /* Trims multiple blanks in line */
       for (j=0; line[j]!='\0';j++){
         line[j]=linetmp[j];
       }
     
   
 int main(int argc, char *argv[])      for (j=maxwav;j>=1;j--){
 {        cutv(stra, strb,line,' '); 
   int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);        if(strb[0]=='.') { /* Missing status */
   int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;          lval=-1;
   int linei, month, year,iout;        }else{
   int jj, ll, li, lj, lk, imk;          errno=0;
   int numlinepar=0; /* Current linenumber of parameter file */          lval=strtol(strb,&endptr,10); 
   int itimes;        /*        if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
   int NDIM=2;          if( strb[0]=='\0' || (*endptr != '\0')){
   int vpopbased=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);
             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);
             return 1;
           }
         }
         s[j][i]=lval;
         
         strcpy(line,stra);
         cutv(stra, strb,line,' ');
         if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
         }
         else  if(iout=sscanf(strb,"%s.") != 0){
           month=99;
           year=9999;
         }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);
           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);
           return 1;
         }
         anint[j][i]= (double) year; 
         mint[j][i]= (double)month; 
         strcpy(line,stra);
       } /* ENd Waves */
       
       cutv(stra, strb,line,' '); 
       if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
       }
       else  if(iout=sscanf(strb,"%s.",dummy) != 0){
         month=99;
         year=9999;
       }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);
           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);
           return 1;
       }
       andc[i]=(double) year; 
       moisdc[i]=(double) month; 
       strcpy(line,stra);
       
       cutv(stra, strb,line,' '); 
       if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
       }
       else  if(iout=sscanf(strb,"%s.") != 0){
         month=99;
         year=9999;
       }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);
         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);
           return 1;
       }
       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);
         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);
           return 1;
   
   char ca[32], cb[32], cc[32];      }
   char dummy[]="                         ";      annais[i]=(double)(year);
   /*  FILE *fichtm; *//* Html File */      moisnais[i]=(double)(month); 
   /* FILE *ficgp;*/ /*Gnuplot File */      strcpy(line,stra);
   struct stat info;      
   double agedeb, agefin,hf;      cutv(stra, strb,line,' '); 
   double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;      errno=0;
       dval=strtod(strb,&endptr); 
       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);
         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);
         fflush(ficlog);
         return 1;
       }
       weight[i]=dval; 
       strcpy(line,stra);
       
       for (j=ncovcol;j>=1;j--){
         cutv(stra, strb,line,' '); 
         if(strb[0]=='.') { /* Missing status */
           lval=-1;
         }else{
           errno=0;
           lval=strtol(strb,&endptr,10); 
           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);
             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);
             return 1;
           }
         }
         if(lval <-1 || lval >1){
           printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\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 example, for multinomial values like 1, 2 and 3,\n \
    build V1=0 V2=0 for the reference value (1),\n \
           V1=1 V2=0 for (2) \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 \
    Exiting.\n",lval,linei, i,line,j);
           fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\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 example, for multinomial values like 1, 2 and 3,\n \
    build V1=0 V2=0 for the reference value (1),\n \
           V1=1 V2=0 for (2) \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 \
    Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
           return 1;
         }
         covar[j][i]=(double)(lval);
         strcpy(line,stra);
       }  
       lstra=strlen(stra);
        
       if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */
         stratrunc = &(stra[lstra-9]);
         num[i]=atol(stratrunc);
       }
       else
         num[i]=atol(stra);
       /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
         printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
       
       i=i+1;
     } /* End loop reading  data */
   
   double fret;    *imax=i-1; /* Number of individuals */
   double **xi,tmp,delta;    fclose(fic);
    
     return (0);
     endread:
       printf("Exiting readdata: ");
       fclose(fic);
       return (1);
   
   double dum; /* Dummy variable */  
   double ***p3mat;  
   double ***mobaverage;  
   int *indx;  
   char line[MAXLINE], linepar[MAXLINE];  
   char linetmp[MAXLINE];  
     char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];  
   char pathr[MAXLINE], pathimach[MAXLINE];   
   char **bp, *tok, *val; /* pathtot */  
   int firstobs=1, lastobs=10;  
   int sdeb, sfin; /* Status at beginning and end */  
   int c,  h , cpt,l;  
   int ju,jl, mi;  
   int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;  
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,*tab;   
   int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */  
   int mobilav=0,popforecast=0;  
   int hstepm, nhstepm;  
   int agemortsup;  
   float  sumlpop=0.;  
   double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000;  
   double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000;  
   
   double bage, fage, age, agelim, agebase;  
   double ftolpl=FTOL;  
   double **prlim;  
   double *severity;  
   double ***param; /* Matrix of parameters */  
   double  *p;  
   double **matcov; /* Matrix of covariance */  
   double ***delti3; /* Scale */  
   double *delti; /* Scale */  
   double ***eij, ***vareij;  
   double **varpl; /* Variances of prevalence limits by age */  
   double *epj, vepp;  
   double kk1, kk2;  
   double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;  
   double **ximort;  
   char *alph[]={"a","a","b","c","d","e"}, str[4];  
   int *dcwave;  
   
   char z[1]="c", occ;  }
   
   char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];  int decodemodel ( char model[], int lastobs)
   char  *strt, strtend[80];  {
   char *stratrunc;    int i, j, k;
   int lstra;    int i1, j1, k1, k2;
     char modelsav[80];
      char stra[80], strb[80], strc[80], strd[80],stre[80];
   
   long total_usecs;    if (strlen(model) >1){ /* If there is at least 1 covariate */
        j=0, j1=0, k1=1, k2=1;
 /*   setlocale (LC_ALL, ""); */      j=nbocc(model,'+'); /* j=Number of '+' */
 /*   bindtextdomain (PACKAGE, LOCALEDIR); */      j1=nbocc(model,'*'); /* j1=Number of '*' */
 /*   textdomain (PACKAGE); */      cptcovn=j+1; /* Number of covariates V1+V2*age+V3 =>(2 plus signs) + 1=3 
                     but the covariates which are product must be computed and stored. */
       cptcovprod=j1; /*Number of products  V1*V2 +v3*age = 2 */
       
       strcpy(modelsav,model); 
       if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){
         printf("Error. Non available option model=%s ",model);
         fprintf(ficlog,"Error. Non available option model=%s ",model);fflush(ficlog);
         return 1;
       }
       
       /* 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 */
       /*    modelsav=V3*age+V2+V1+V4 strb=V3*age stra=V2+V1+V4 
           i=1 Tvar[1]=3 Tage[1]=1  
           i=2 Tvar[2]=2
           i=3 Tvar[3]=1
           i=4 Tvar[4]= 4
           i=5 Tvar[5]
         for (k=1; k<=cptcovn;k++) {
           cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
        */
       for(k=1; k<=(j+1);k++){
         cutv(strb,stra,modelsav,'+'); /* keeps in strb after the first '+' 
                                        modelsav=V3*age+V2+V1+V4 strb=V3*age stra=V2+V1+V4 
                                       */ 
         /* if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav);*/ /* and analyzes it */
         /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
         /*scanf("%d",i);*/
         if (strchr(strb,'*')) {  /* Model includes a product V3*age+V2+V1+V4 strb=V3*age */
           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 */
             cptcovprod--;
             cutv(strb,stre,strd,'V'); /* stre="V3" */
             Tvar[k]=atoi(stre);  /* V1+V3*age+V2 Tvar[2]=3, and Tvar[3]=2 */
             cptcovage++; /* Sums the number of covariates which include age as a product */
             Tage[cptcovage]=k;  /* Tage[1] =2 */
             /*printf("stre=%s ", stre);*/
           }
           else if (strcmp(strd,"age")==0) { /* or age*Vn */
             cptcovprod--;
             cutv(strb,stre,strc,'V');
             Tvar[k]=atoi(stre);
             cptcovage++;
             Tage[cptcovage]=k;
           }
           else {  /* Age is not in the model V1+V3*V2+V2  strb=V3*V2*/
             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. 
                                     If already ncovcol=2 and model=V2*V1 Tvar[1]=2+1 and Tvar[2]=2+2 etc */
             cutv(strb,strc,strd,'V'); /* strd was Vm, strc is m */
             Tprod[k1]=k;  /* Tprod[1]  */
             Tvard[k1][1]=atoi(strc); /* m*/
             Tvard[k1][2]=atoi(stre); /* n */
             Tvar[cptcovn+k2]=Tvard[k1][1];
             Tvar[cptcovn+k2+1]=Tvard[k1][2]; 
             for (i=1; i<=lastobs;i++) /* Computes the new covariate which is a product of covar[n][i]* covar[m][i]
                                        and is stored at ncovol+k1 */
               covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
             k1++;
             k2=k2+2;
           }
         }
         else { /* no more sum */
           /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
          /*  scanf("%d",i);*/
           cutv(strd,strc,strb,'V');
           Tvar[i]=atoi(strc);
         }
         strcpy(modelsav,stra);  /* modelsav=V2+V3*age+V1+V4 strb=V3*age+V1+V4 */ 
         /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
           scanf("%d",i);*/
       } /* end of loop + */
     } /* end model */
     
     /*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products.
       If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/
   
     /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
     printf("cptcovprod=%d ", cptcovprod);
     fprintf(ficlog,"cptcovprod=%d ", cptcovprod);
   
     scanf("%d ",i);*/
   
   
     return (0);
     endread:
       printf("Exiting decodemodel: ");
       return (1);
   }
   
   calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )
   {
     int i, m;
   
     for (i=1; i<=imx; i++) {
       for(m=2; (m<= maxwav); m++) {
         if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){
           anint[m][i]=9999;
           s[m][i]=-1;
         }
         if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){
           *nberr++;
           printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);
           fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);
           s[m][i]=-1;
         }
         if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
           *nberr++;
           printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); 
           fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); 
           s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */
         }
       }
     }
   
     for (i=1; i<=imx; i++)  {
       agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
       for(m=firstpass; (m<= lastpass); m++){
         if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){
           if (s[m][i] >= nlstate+1) {
             if(agedc[i]>0)
               if((int)moisdc[i]!=99 && (int)andc[i]!=9999)
                 agev[m][i]=agedc[i];
             /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
               else {
                 if ((int)andc[i]!=9999){
                   nbwarn++;
                   printf("Warning negative age at death: %ld line:%d\n",num[i],i);
                   fprintf(ficlog,"Warning negative age at death: %ld line:%d\n",num[i],i);
                   agev[m][i]=-1;
                 }
               }
           }
           else if(s[m][i] !=9){ /* Standard case, age in fractional
                                    years but with the precision of a month */
             agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
             if((int)mint[m][i]==99 || (int)anint[m][i]==9999)
               agev[m][i]=1;
             else if(agev[m][i] < *agemin){ 
               *agemin=agev[m][i];
               printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], *agemin);
             }
             else if(agev[m][i] >*agemax){
               *agemax=agev[m][i];
               /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/
             }
             /*agev[m][i]=anint[m][i]-annais[i];*/
             /*     agev[m][i] = age[i]+2*m;*/
           }
           else { /* =9 */
             agev[m][i]=1;
             s[m][i]=-1;
           }
         }
         else /*= 0 Unknown */
           agev[m][i]=1;
       }
       
     }
     for (i=1; i<=imx; i++)  {
       for(m=firstpass; (m<=lastpass); m++){
         if (s[m][i] > (nlstate+ndeath)) {
           *nberr++;
           printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);     
           fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);     
           return 1;
         }
       }
     }
   
     /*for (i=1; i<=imx; i++){
     for (m=firstpass; (m<lastpass); m++){
        printf("%ld %d %.lf %d %d\n", num[i],(covar[1][i]),agev[m][i],s[m][i],s[m+1][i]);
   }
   
   }*/
   
   
     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); 
   
     return (0);
     endread:
       printf("Exiting calandcheckages: ");
       return (1);
   }
   
   
   /***********************************************/
   /**************** Main Program *****************/
   /***********************************************/
   
   int main(int argc, char *argv[])
   {
   #ifdef GSL
     const gsl_multimin_fminimizer_type *T;
     size_t iteri = 0, it;
     int rval = GSL_CONTINUE;
     int status = GSL_SUCCESS;
     double ssval;
   #endif
     int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
     int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;
     int linei, month, year,iout;
     int jj, ll, li, lj, lk, imk;
     int numlinepar=0; /* Current linenumber of parameter file */
     int itimes;
     int NDIM=2;
     int vpopbased=0;
   
     char ca[32], cb[32], cc[32];
     /*  FILE *fichtm; *//* Html File */
     /* FILE *ficgp;*/ /*Gnuplot File */
     struct stat info;
     double agedeb, agefin,hf;
     double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
   
     double fret;
     double **xi,tmp,delta;
   
     double dum; /* Dummy variable */
     double ***p3mat;
     double ***mobaverage;
     int *indx;
     char line[MAXLINE], linepar[MAXLINE];
     char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];
     char pathr[MAXLINE], pathimach[MAXLINE]; 
     char **bp, *tok, *val; /* pathtot */
     int firstobs=1, lastobs=10;
     int sdeb, sfin; /* Status at beginning and end */
     int c,  h , cpt,l;
     int ju,jl, mi;
     int i1,j1, jk,aa,bb, stepsize, ij;
     int jnais,jdc,jint4,jint1,jint2,jint3,*tab; 
     int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
     int mobilav=0,popforecast=0;
     int hstepm, nhstepm;
     int agemortsup;
     float  sumlpop=0.;
     double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000;
     double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000;
   
     double bage, fage, age, agelim, agebase;
     double ftolpl=FTOL;
     double **prlim;
     double ***param; /* Matrix of parameters */
     double  *p;
     double **matcov; /* Matrix of covariance */
     double ***delti3; /* Scale */
     double *delti; /* Scale */
     double ***eij, ***vareij;
     double **varpl; /* Variances of prevalence limits by age */
     double *epj, vepp;
     double kk1, kk2;
     double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
     double **ximort;
     char *alph[]={"a","a","b","c","d","e"}, str[4];
     int *dcwave;
   
     char z[1]="c", occ;
   
     /*char  *strt;*/
     char strtend[80];
   
     long total_usecs;
    
   /*   setlocale (LC_ALL, ""); */
   /*   bindtextdomain (PACKAGE, LOCALEDIR); */
   /*   textdomain (PACKAGE); */
 /*   setlocale (LC_CTYPE, ""); */  /*   setlocale (LC_CTYPE, ""); */
 /*   setlocale (LC_MESSAGES, ""); */  /*   setlocale (LC_MESSAGES, ""); */
   
Line 4671  int main(int argc, char *argv[]) Line 5113  int main(int argc, char *argv[])
   
         
   covar=matrix(0,NCOVMAX,1,n);     covar=matrix(0,NCOVMAX,1,n); 
   cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/    cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/
   if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;    /* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5
   /* where is ncovprod ?*/       v1+v2*age+v2*v3 makes cptcovn = 3
   ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/    */
     if (strlen(model)>1) 
       cptcovn=nbocc(model,'+')+1;
     /* ncovprod */
     ncovmodel=2+cptcovn; /*Number of variables including intercept and age = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/
   nvar=ncovmodel-1; /* Suppressing age as a basic covariate */    nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
   nforce= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */    nforce= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */
   npar= nforce*ncovmodel; /* Number of parameters like aij*/    npar= nforce*ncovmodel; /* Number of parameters like aij*/
Line 4856  run imach with mle=-1 to get a correct t Line 5302  run imach with mle=-1 to get a correct t
     fprintf(ficres,"#%s\n",version);      fprintf(ficres,"#%s\n",version);
   }    /* End of mle != -3 */    }    /* End of mle != -3 */
   
   /*-------- data file ----------*/  
   if((fic=fopen(datafile,"r"))==NULL)    {  
     printf("Problem while opening datafile: %s\n", datafile);goto end;  
     fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);goto end;  
   }  
   
   n= lastobs;    n= lastobs;
   severity = vector(1,maxwav);  
   outcome=imatrix(1,maxwav+1,1,n);  
   num=lvector(1,n);    num=lvector(1,n);
   moisnais=vector(1,n);    moisnais=vector(1,n);
   annais=vector(1,n);    annais=vector(1,n);
Line 4880  run imach with mle=-1 to get a correct t Line 5319  run imach with mle=-1 to get a correct t
   tab=ivector(1,NCOVMAX);    tab=ivector(1,NCOVMAX);
   ncodemax=ivector(1,8);    ncodemax=ivector(1,8);
   
   i=1;    /* Reads data from file datafile */
   linei=0;    if (readdata(datafile, firstobs, lastobs, &imx)==1)
   while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) {      goto end;
     linei=linei+1;  
     for(j=strlen(line); j>=0;j--){  /* Untabifies line */  
       if(line[j] == '\t')  
         line[j] = ' ';  
     }  
     for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10)||(line[j]==13);j--){  
       ;  
     };  
     line[j+1]=0;  /* Trims blanks at end of line */  
     if(line[0]=='#'){  
       fprintf(ficlog,"Comment line\n%s\n",line);  
       printf("Comment line\n%s\n",line);  
       continue;  
     }  
     trimbb(linetmp,line); /* Trims multiple blanks in line */  
     for (j=0; line[j]!='\0';j++){  
       line[j]=linetmp[j];  
     }  
     
   
     for (j=maxwav;j>=1;j--){  
       cutv(stra, strb,line,' ');   
       if(strb[0]=='.') { /* Missing status */  
         lval=-1;  
       }else{  
         errno=0;  
         lval=strtol(strb,&endptr,10);   
       /*        if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/  
         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);  
           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);  
           goto end;  
         }  
       }  
       s[j][i]=lval;  
         
       strcpy(line,stra);  
       cutv(stra, strb,line,' ');  
       if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){  
       }  
       else  if(iout=sscanf(strb,"%s.") != 0){  
         month=99;  
         year=9999;  
       }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);  
         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);  
         goto end;  
       }  
       anint[j][i]= (double) year;   
       mint[j][i]= (double)month;   
       strcpy(line,stra);  
     } /* ENd Waves */  
       
     cutv(stra, strb,line,' ');   
     if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){  
     }  
     else  if(iout=sscanf(strb,"%s.",dummy) != 0){  
       month=99;  
       year=9999;  
     }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);  
         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);  
         goto end;  
     }  
     andc[i]=(double) year;   
     moisdc[i]=(double) month;   
     strcpy(line,stra);  
       
     cutv(stra, strb,line,' ');   
     if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){  
     }  
     else  if(iout=sscanf(strb,"%s.") != 0){  
       month=99;  
       year=9999;  
     }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);  
       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);  
         goto end;  
     }  
     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);  
       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);  
         goto end;  
   
     }  
     annais[i]=(double)(year);  
     moisnais[i]=(double)(month);   
     strcpy(line,stra);  
       
     cutv(stra, strb,line,' ');   
     errno=0;  
     dval=strtod(strb,&endptr);   
     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);  
       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);  
       fflush(ficlog);  
       goto end;  
     }  
     weight[i]=dval;   
     strcpy(line,stra);  
       
     for (j=ncovcol;j>=1;j--){  
       cutv(stra, strb,line,' ');   
       if(strb[0]=='.') { /* Missing status */  
         lval=-1;  
       }else{  
         errno=0;  
         lval=strtol(strb,&endptr,10);   
         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);  
           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);  
           goto end;  
         }  
       }  
       if(lval <-1 || lval >1){  
         printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\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 example, for multinomial values like 1, 2 and 3,\n \  
  build V1=0 V2=0 for the reference value (1),\n \  
         V1=1 V2=0 for (2) \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 \  
  Exiting.\n",lval,linei, i,line,j);  
         fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\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 example, for multinomial values like 1, 2 and 3,\n \  
  build V1=0 V2=0 for the reference value (1),\n \  
         V1=1 V2=0 for (2) \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 \  
  Exiting.\n",lval,linei, i,line,j);fflush(ficlog);  
         goto end;  
       }  
       covar[j][i]=(double)(lval);  
       strcpy(line,stra);  
     }    
     lstra=strlen(stra);  
        
     if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */  
       stratrunc = &(stra[lstra-9]);  
       num[i]=atol(stratrunc);  
     }  
     else  
       num[i]=atol(stra);  
     /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){  
       printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/  
       
     i=i+1;  
   } /* End loop reading  data */  
   fclose(fic);  
   /* printf("ii=%d", ij);  
      scanf("%d",i);*/  
   imx=i-1; /* Number of individuals */  
   
   /* for (i=1; i<=imx; i++){  
     if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;  
     if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;  
     if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;  
     }*/  
    /*  for (i=1; i<=imx; i++){  
      if (s[4][i]==9)  s[4][i]=-1;   
      printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));}*/  
     
   /* for (i=1; i<=imx; i++) */  
    
    /*if ((s[3][i]==3) ||  (s[4][i]==3)) weight[i]=0.08;  
      else weight[i]=1;*/  
   
   /* 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 */    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 */
Line 5059  run imach with mle=-1 to get a correct t Line 5329  run imach with mle=-1 to get a correct t
   Tvaraff=ivector(1,15);     Tvaraff=ivector(1,15); 
   Tvard=imatrix(1,15,1,2);    Tvard=imatrix(1,15,1,2);
   Tage=ivector(1,15);          Tage=ivector(1,15);      
      
   if (strlen(model) >1){ /* If there is at least 1 covariate */  
     j=0, j1=0, k1=1, k2=1;  
     j=nbocc(model,'+'); /* j=Number of '+' */  
     j1=nbocc(model,'*'); /* j1=Number of '*' */  
     cptcovn=j+1; /* Number of covariates V1+V2+V3 =>2+1=3 */  
     cptcovprod=j1; /*Number of products  V1*V2 =1 */  
       
     strcpy(modelsav,model);   
     if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){  
       printf("Error. Non available option model=%s ",model);  
       fprintf(ficlog,"Error. Non available option model=%s ",model);fflush(ficlog);  
       goto end;  
     }  
       
     /* 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 */  
     for(i=(j+1); i>=1;i--){  
       cutv(stra,strb,modelsav,'+'); /* keeps in strb after the first '+'   
                                      modelsav=V2+V3*age+V1+V4 strb=V3*age+V1+V4   
                                      stra=V2  
                                     */   
       if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */  
       /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/  
       /*scanf("%d",i);*/  
       if (strchr(strb,'*')) {  /* Model includes a product V1+V3*age+V2 strb=V3*age*/  
         cutv(strd,strc,strb,'*'); /* strd*strc  Vm*Vn: V3*age strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */  
         if (strcmp(strc,"age")==0) { /* Vn*age */  
           cptcovprod--;  
           cutv(strb,stre,strd,'V');  
           Tvar[i]=atoi(stre);  /* V1+V3*age+V2 Tvar[2]=3 */  
           cptcovage++; /* Sums the number of covariates including age as a product */  
           Tage[cptcovage]=i;  /* Tage[1] =2 */  
           /*printf("stre=%s ", stre);*/  
         }  
         else if (strcmp(strd,"age")==0) { /* or age*Vn */  
           cptcovprod--;  
           cutv(strb,stre,strc,'V');  
           Tvar[i]=atoi(stre);  
           cptcovage++;  
           Tage[cptcovage]=i;  
         }  
         else {  /* Age is not in the model V1+V3*V2+V2  strb=V3*V2*/  
           cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/  
           Tvar[i]=ncovcol+k1;  /* find 'n' in Vn and stores in Tvar.   
                                   If already ncovcol=2 and model=V2*V1 Tvar[1]=2+1 and Tvar[2]=2+2 etc */  
           cutv(strb,strc,strd,'V'); /* strd was Vm, strc is m */  
           Tprod[k1]=i;  /* Tprod[1]  */  
           Tvard[k1][1]=atoi(strc); /* m*/  
           Tvard[k1][2]=atoi(stre); /* n */  
           Tvar[cptcovn+k2]=Tvard[k1][1];  
           Tvar[cptcovn+k2+1]=Tvard[k1][2];   
           for (k=1; k<=lastobs;k++)   
             covar[ncovcol+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];  
           k1++;  
           k2=k2+2;  
         }  
       }  
       else { /* no more sum */  
         /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/  
        /*  scanf("%d",i);*/  
       cutv(strd,strc,strb,'V');  
       Tvar[i]=atoi(strc);  
       }  
       strcpy(modelsav,stra);  /* modelsav=V2+V3*age+V1+V4 strb=V3*age+V1+V4 */   
       /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);  
         scanf("%d",i);*/  
     } /* end of loop + */  
   } /* end model */  
     
   /*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products.  
     If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/  
   
   /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);    if(decodemodel(model, lastobs) == 1)
   printf("cptcovprod=%d ", cptcovprod);      goto end;
   fprintf(ficlog,"cptcovprod=%d ", cptcovprod);  
   
   scanf("%d ",i);*/  
   
     /*  if(mle==1){*/      /*  if(mle==1){*/
   if (weightopt != 1) { /* Maximisation without weights*/    if (weightopt != 1) { /* Maximisation without weights*/
     for(i=1;i<=n;i++) weight[i]=1.0;      for(i=1;i<=n;i++) weight[i]=1.0;
   }    }
   
     /*-calculation of age at interview from date of interview and age at death -*/      /*-calculation of age at interview from date of interview and age at death -*/
   agev=matrix(1,maxwav,1,imx);    agev=matrix(1,maxwav,1,imx);
   
   for (i=1; i<=imx; i++) {    if(calandcheckages(imx, maxwav, &agemin, &agemax, &nberr, &nbwarn) == 1)
     for(m=2; (m<= maxwav); m++) {      goto end;
       if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){  
         anint[m][i]=9999;  
         s[m][i]=-1;  
       }  
       if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){  
         nberr++;  
         printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);  
         fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);  
         s[m][i]=-1;  
       }  
       if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){  
         nberr++;  
         printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]);   
         fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]);   
         s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */  
       }  
     }  
   }  
   
   for (i=1; i<=imx; i++)  {  
     agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);  
     for(m=firstpass; (m<= lastpass); m++){  
       if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){  
         if (s[m][i] >= nlstate+1) {  
           if(agedc[i]>0)  
             if((int)moisdc[i]!=99 && (int)andc[i]!=9999)  
               agev[m][i]=agedc[i];  
           /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/  
             else {  
               if ((int)andc[i]!=9999){  
                 nbwarn++;  
                 printf("Warning negative age at death: %ld line:%d\n",num[i],i);  
                 fprintf(ficlog,"Warning negative age at death: %ld line:%d\n",num[i],i);  
                 agev[m][i]=-1;  
               }  
             }  
         }  
         else if(s[m][i] !=9){ /* Standard case, age in fractional  
                                  years but with the precision of a month */  
           agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);  
           if((int)mint[m][i]==99 || (int)anint[m][i]==9999)  
             agev[m][i]=1;  
           else if(agev[m][i] <agemin){   
             agemin=agev[m][i];  
             printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);  
           }  
           else if(agev[m][i] >agemax){  
             agemax=agev[m][i];  
             /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/  
           }  
           /*agev[m][i]=anint[m][i]-annais[i];*/  
           /*     agev[m][i] = age[i]+2*m;*/  
         }  
         else { /* =9 */  
           agev[m][i]=1;  
           s[m][i]=-1;  
         }  
       }  
       else /*= 0 Unknown */  
         agev[m][i]=1;  
     }  
       
   }  
   for (i=1; i<=imx; i++)  {  
     for(m=firstpass; (m<=lastpass); m++){  
       if (s[m][i] > (nlstate+ndeath)) {  
         nberr++;  
         printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);       
         fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);       
         goto end;  
       }  
     }  
   }  
   
   /*for (i=1; i<=imx; i++){  
   for (m=firstpass; (m<lastpass); m++){  
      printf("%ld %d %.lf %d %d\n", num[i],(covar[1][i]),agev[m][i],s[m][i],s[m+1][i]);  
 }  
   
 }*/  
   
   
   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);   
   
   agegomp=(int)agemin;    agegomp=(int)agemin;
   free_vector(severity,1,maxwav);  
   free_imatrix(outcome,1,maxwav+1,1,n);  
   free_vector(moisnais,1,n);    free_vector(moisnais,1,n);
   free_vector(annais,1,n);    free_vector(annais,1,n);
   /* free_matrix(mint,1,maxwav,1,n);    /* free_matrix(mint,1,maxwav,1,n);
Line 5267  run imach with mle=-1 to get a correct t Line 5378  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;codtab[h][k]=j;codtab[h][Tvar[k]]=j;              h=1;
               codtab[h][k]=j;
               codtab[h][Tvar[k]]=j;
             }
           printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d 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 5366  Interval (in months) between two waves: Line 5480  Interval (in months) between two waves:
   globpr=0; /* To get the number ipmx of contributions and the sum of weights*/    globpr=0; /* To get the number ipmx of contributions and the sum of weights*/
   
   if (mle==-3){    if (mle==-3){
     ximort=matrix(1,NDIM,1,NDIM);      ximort=matrix(1,NDIM,1,NDIM); 
   /*     ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */
     cens=ivector(1,n);      cens=ivector(1,n);
     ageexmed=vector(1,n);      ageexmed=vector(1,n);
     agecens=vector(1,n);      agecens=vector(1,n);
Line 5408  Interval (in months) between two waves: Line 5523  Interval (in months) between two waves:
     /*printf("%lf %lf", p[1], p[2]);*/      /*printf("%lf %lf", p[1], p[2]);*/
           
           
   #ifdef GSL
       printf("GSL optimization\n");  fprintf(ficlog,"Powell\n");
   #elsedef
     printf("Powell\n");  fprintf(ficlog,"Powell\n");      printf("Powell\n");  fprintf(ficlog,"Powell\n");
   #endif
     strcpy(filerespow,"pow-mort");       strcpy(filerespow,"pow-mort"); 
     strcat(filerespow,fileres);      strcat(filerespow,fileres);
     if((ficrespow=fopen(filerespow,"w"))==NULL) {      if((ficrespow=fopen(filerespow,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", filerespow);        printf("Problem with resultfile: %s\n", filerespow);
       fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);        fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
     }      }
   #ifdef GSL
       fprintf(ficrespow,"# GSL optimization\n# iter -2*LL");
   #elsedef
     fprintf(ficrespow,"# Powell\n# iter -2*LL");      fprintf(ficrespow,"# Powell\n# iter -2*LL");
   #endif
     /*  for (i=1;i<=nlstate;i++)      /*  for (i=1;i<=nlstate;i++)
         for(j=1;j<=nlstate+ndeath;j++)          for(j=1;j<=nlstate+ndeath;j++)
         if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);          if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
     */      */
     fprintf(ficrespow,"\n");      fprintf(ficrespow,"\n");
   #ifdef GSL
       /* gsl starts here */ 
       T = gsl_multimin_fminimizer_nmsimplex;
       gsl_multimin_fminimizer *sfm = NULL;
       gsl_vector *ss, *x;
       gsl_multimin_function minex_func;
   
       /* Initial vertex size vector */
       ss = gsl_vector_alloc (NDIM);
       
       if (ss == NULL){
         GSL_ERROR_VAL ("failed to allocate space for ss", GSL_ENOMEM, 0);
       }
       /* Set all step sizes to 1 */
       gsl_vector_set_all (ss, 0.001);
   
       /* Starting point */
       
       x = gsl_vector_alloc (NDIM);
       
       if (x == NULL){
         gsl_vector_free(ss);
         GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
       }
     
       /* Initialize method and iterate */
       /*     p[1]=0.0268; p[NDIM]=0.083; */
   /*     gsl_vector_set(x, 0, 0.0268); */
   /*     gsl_vector_set(x, 1, 0.083); */
       gsl_vector_set(x, 0, p[1]);
       gsl_vector_set(x, 1, p[2]);
   
       minex_func.f = &gompertz_f;
       minex_func.n = NDIM;
       minex_func.params = (void *)&p; /* ??? */
           
     powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);      sfm = gsl_multimin_fminimizer_alloc (T, NDIM);
       gsl_multimin_fminimizer_set (sfm, &minex_func, x, ss);
       
       printf("Iterations beginning .....\n\n");
       printf("Iter. #    Intercept       Slope     -Log Likelihood     Simplex size\n");
   
       iteri=0;
       while (rval == GSL_CONTINUE){
         iteri++;
         status = gsl_multimin_fminimizer_iterate(sfm);
         
         if (status) printf("error: %s\n", gsl_strerror (status));
         fflush(0);
         
         if (status) 
           break;
         
         rval = gsl_multimin_test_size (gsl_multimin_fminimizer_size (sfm), 1e-6);
         ssval = gsl_multimin_fminimizer_size (sfm);
         
         if (rval == GSL_SUCCESS)
           printf ("converged to a local maximum at\n");
         
         printf("%5d ", iteri);
         for (it = 0; it < NDIM; it++){
           printf ("%10.5f ", gsl_vector_get (sfm->x, it));
         }
         printf("f() = %-10.5f ssize = %.7f\n", sfm->fval, ssval);
       }
       
       printf("\n\n Please note: Program should be run many times with varying starting points to detemine global maximum\n\n");
       
       gsl_vector_free(x); /* initial values */
       gsl_vector_free(ss); /* inital step size */
       for (it=0; it<NDIM; it++){
         p[it+1]=gsl_vector_get(sfm->x,it);
         fprintf(ficrespow," %.12lf", p[it]);
       }
       gsl_multimin_fminimizer_free (sfm); /* p *(sfm.x.data) et p *(sfm.x.data+1)  */
   #endif
   #ifdef POWELL
        powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);
   #endif  
     fclose(ficrespow);      fclose(ficrespow);
           
     hesscov(matcov, p, NDIM, delti, 1e-4, gompertz);       hesscov(matcov, p, NDIM, delti, 1e-4, gompertz); 
Line 5483  Interval (in months) between two waves: Line 5683  Interval (in months) between two waves:
     free_vector(lsurv,1,AGESUP);      free_vector(lsurv,1,AGESUP);
     free_vector(lpop,1,AGESUP);      free_vector(lpop,1,AGESUP);
     free_vector(tpop,1,AGESUP);      free_vector(tpop,1,AGESUP);
   #ifdef GSL
       free_ivector(cens,1,n);
       free_vector(agecens,1,n);
       free_ivector(dcwave,1,n);
       free_matrix(ximort,1,NDIM,1,NDIM);
   #endif
   } /* Endof if mle==-3 */    } /* Endof if mle==-3 */
       
   else{ /* For mle >=1 */    else{ /* For mle >=1 */
Removed from v.1.135  
changed lines
  Added in v.1.136

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