--- imach/src/imach.c 2009/05/26 06:44:34 1.130 +++ imach/src/imach.c 2009/06/20 16:22:47 1.131 @@ -1,6 +1,9 @@ -/* $Id: imach.c,v 1.130 2009/05/26 06:44:34 brouard Exp $ +/* $Id: imach.c,v 1.131 2009/06/20 16:22:47 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.131 2009/06/20 16:22:47 brouard + Some dimensions resccaled + Revision 1.130 2009/05/26 06:44:34 brouard (Module): Max Covariate is now set to 20 instead of 8. A lot of cleaning with variables initialized to 0. Trying to make @@ -372,7 +375,7 @@ extern int errno; #define GLOCK_ERROR_NOPATH -1 /* empty path */ #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ -#define MAXPARM 30 /* Maximum number of parameters for the optimization */ +#define MAXPARM 128 /* Maximum number of parameters for the optimization */ #define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */ #define NINTERVMAX 8 @@ -394,11 +397,11 @@ extern int errno; #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.130 2009/05/26 06:44:34 brouard Exp $ */ +/* $Id: imach.c,v 1.131 2009/06/20 16:22:47 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.98i, June 2006, INED-EUROREVES-Institut de longevite "; -char fullversion[]="$Revision: 1.130 $ $Date: 2009/05/26 06:44:34 $"; +char version[]="Imach version 0.98k, June 2006, INED-EUROREVES-Institut de longevite "; +char fullversion[]="$Revision: 1.131 $ $Date: 2009/06/20 16:22:47 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -521,7 +524,7 @@ double dateintmean=0; double *weight; int **s; /* Status */ double *agedc, **covar, idx; -int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff; +int **nbcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff; double *lsurv, *lpop, *tpop; double ftol=FTOL; /* Tolerance for computing Max Likelihood */ @@ -1171,7 +1174,7 @@ double **prevalim(double **prlim, int nl int i, ii,j,k; double min, max, maxmin, maxmax,sumnew=0.; double **matprod2(); - double **out, cov[NCOVMAX], **pmij(); + double **out, cov[NCOVMAX+1], **pmij(); double **newm; double agefin, delaymax=50 ; /* Max number of years to converge */ @@ -1253,10 +1256,14 @@ double **pmij(double **ps, double *cov, for(i=1; i<= nlstate; i++){ s1=0; - for(j=1; j cptcode) cptcode=ij; /* getting the maximum of covariable + if (ij > cptcode) cptcode=ij; /* getting the maximum value of the modality of the covariate (should be 0 or 1 now) Tvar[j]. If V=sex and male is 0 and female is 1, then cptcode=1.*/ } - for (i=0; i<=cptcode; i++) { - if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariates. In fact ncodemax[j]=2 (dichotom. variables) but it can be more */ - } + for (i=0; i<=cptcode; i++) { /* i=-1 ?*/ + if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariate. In fact ncodemax[j]=2 (dichotom. variables only) but it can be more */ + } /* Ndum[-1] number of undefined modalities */ ij=1; - for (i=1; i<=ncodemax[j]; i++) { - for (k=0; k<= maxncov; k++) { - if (Ndum[k] != 0) { - nbcode[Tvar[j]][ij]=k; - /* store the modality in an array. k is a modality. If we have model=V1+V1*sex then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ - + for (i=1; i<=ncodemax[j]; i++) { /* i= 1 to 2 */ + for (k=0; k<= maxncov; k++) { /* k=-1 ?*/ + 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. + k is a modality. If we have model=V1+V1*sex + then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ ij++; } if (ij > ncodemax[j]) break; @@ -2475,9 +2483,9 @@ void tricode(int *Tvar, int **nbcode, in for (k=0; k< maxncov; k++) Ndum[k]=0; - for (i=1; i<=ncovmodel-2; i++) { + 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.*/ - ij=Tvar[i]; + ij=Tvar[i]; /* Tvar might be -1 if status was unknown */ Ndum[ij]++; } @@ -2488,8 +2496,8 @@ void tricode(int *Tvar, int **nbcode, in ij++; } } - - cptcoveff=ij-1; /*Number of simple covariates*/ + ij--; + cptcoveff=ij; /*Number of simple covariates*/ } /*********** Health Expectancies ****************/ @@ -3087,9 +3095,9 @@ void varevsij(char optionfilefiname[], d free_vector(gmp,nlstate+1,nlstate+ndeath); free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath); free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/ - fprintf(ficgp,"\nset noparametric;set nolabel; set ter png small;set size 0.65, 0.65"); + fprintf(ficgp,"\nunset parametric;unset label; set ter png small;set size 0.65, 0.65"); /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */ - fprintf(ficgp,"\n set log y; set nolog x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); + fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); /* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ @@ -3266,7 +3274,7 @@ void varprob(char optionfilefiname[], do fprintf(ficresprobcov,"\n"); fprintf(ficresprobcor,"\n"); */ - xp=vector(1,npar); + xp=vector(1,npar); dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage); @@ -3419,7 +3427,7 @@ To be simple, these graphs help to under /* Confidence intervalle of pij */ /* - fprintf(ficgp,"\nset noparametric;unset label"); + fprintf(ficgp,"\nunset parametric;unset label"); fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\""); fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65"); fprintf(fichtm,"\n
Probability with confidence intervals expressed in year-1 :pijgr%s.png, ",optionfilefiname,optionfilefiname); @@ -3681,7 +3689,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); - else fprintf(ficgp," \%%*lf (\%%*lf)"); + else fprintf(ficgp," \%%*lf (\%%*lf)"); } fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); for (i=1; i<= nlstate ; i ++) { @@ -4631,7 +4639,13 @@ int main(int argc, char *argv[]) ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ - + if(npar >MAXPARM || nlstate >NLSTATEMAX || ndeath >NDEATHMAX || ncovmodel>NCOVMAX){ + printf("Too complex model for current IMaCh: npar=(nlstate+ndeath-1)*nlstate*ncovmodel=%d >= %d(MAXPARM) or nlstate=%d >= %d(NLSTATEMAX) or ndeath=%d >= %d(NDEATHMAX) or ncovmodel=(k+age+#of+signs)=%d(NCOVMAX) >= %d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX); + fprintf(ficlog,"Too complex model for current IMaCh: %d >=%d(MAXPARM) or %d >=%d(NLSTATEMAX) or %d >=%d(NDEATHMAX) or %d(NCOVMAX) >=%d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX); + fflush(stdout); + fclose (ficlog); + goto end; + } delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); delti=delti3[1][1]; /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ @@ -4757,6 +4771,9 @@ run imach with mle=-1 to get a correct t ungetc(c,ficpar); matcov=matrix(1,npar,1,npar); + for(i=1; i <=npar; i++) + for(j=1; j <=npar; j++) matcov[i][j]=0.; + for(i=1; i <=npar; i++){ fscanf(ficpar,"%s",&str); if(mle==1) @@ -4820,7 +4837,7 @@ run imach with mle=-1 to get a correct t for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */ mint=matrix(1,maxwav,1,n); anint=matrix(1,maxwav,1,n); - s=imatrix(1,maxwav+1,1,n); + s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */ tab=ivector(1,NCOVMAX); ncodemax=ivector(1,8); @@ -4844,12 +4861,17 @@ run imach with mle=-1 to get a correct t for (j=maxwav;j>=1;j--){ cutv(stra, strb,line,' '); - errno=0; - lval=strtol(strb,&endptr,10); + 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 '%d' at line number %d %s 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); - exit(1); + 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; @@ -4861,8 +4883,9 @@ run imach with mle=-1 to get a correct t month=99; year=9999; }else{ - printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j); - exit(1); + 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; @@ -4876,8 +4899,9 @@ run imach with mle=-1 to get a correct t month=99; year=9999; }else{ - printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line); - exit(1); + 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; @@ -4891,7 +4915,8 @@ run imach with mle=-1 to get a correct t year=9999; }else{ printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j); - exit(1); + fprintf(ficlog,"Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j);fflush(ficlog); + goto end; } annais[i]=(double)(year); moisnais[i]=(double)(month); @@ -4902,18 +4927,25 @@ run imach with mle=-1 to get a correct t 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); - exit(1); + 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,' '); - errno=0; - lval=strtol(strb,&endptr,10); - if( strb[0]=='\0' || (*endptr != '\0')){ - printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a covar (meaning 0 for the reference or 1). Exiting.\n",lval, linei,i, line); - exit(1); + 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 \ @@ -4925,6 +4957,15 @@ run imach with mle=-1 to get a correct t 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); @@ -4963,7 +5004,7 @@ run imach with mle=-1 to get a correct t else weight[i]=1;*/ /* Calculation of the number of parameters from char model */ - Tvar=ivector(1,15); /* 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 */ Tprod=ivector(1,15); Tvaraff=ivector(1,15); Tvard=imatrix(1,15,1,2); @@ -4979,7 +5020,7 @@ run imach with mle=-1 to get a correct t 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); + fprintf(ficlog,"Error. Non available option model=%s ",model);fflush(ficlog); goto end; } @@ -5162,7 +5203,6 @@ run imach with mle=-1 to get a correct t /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ - Tcode=ivector(1,100); nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); ncodemax[1]=1; if (cptcovn > 0) tricode(Tvar,nbcode,imx); @@ -5172,12 +5212,13 @@ run imach with mle=-1 to get a correct t h=0; m=pow(2,cptcoveff); - for(k=1;k<=cptcoveff; k++){ - for(i=1; i <=(m/pow(2,k));i++){ - for(j=1; j <= ncodemax[k]; j++){ - for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){ + for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */ + for(i=1; i <=(m/pow(2,k));i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */ + for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate */ + 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++; - if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j; + if (h>m) + h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j; printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]); } } @@ -5187,7 +5228,7 @@ run imach with mle=-1 to get a correct t codtab[1][2]=1;codtab[2][2]=2; */ /* for(i=1; i <=m ;i++){ for(k=1; k <=cptcovn; k++){ - printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); + printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); } printf("\n"); } @@ -5215,7 +5256,8 @@ run imach with mle=-1 to get a correct t strcat(optionfilehtm,"-mort"); strcat(optionfilehtm,".htm"); if((fichtm=fopen(optionfilehtm,"w"))==NULL) { - printf("Problem with %s \n",optionfilehtm), exit(0); + printf("Problem with %s \n",optionfilehtm); + exit(0); } strcpy(optionfilehtmcov,optionfilefiname); /* Only for matrix of covariance */ @@ -5394,8 +5436,8 @@ Interval (in months) between two waves: } /* Endof if mle==-3 */ else{ /* For mle >=1 */ - - likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ + globpr=1;/* debug */ + /* likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone);*/ /* Prints the contributions to the likelihood */ printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); for (k=1; k<=npar;k++) printf(" %d %8.5f",k,p[k]); @@ -5667,7 +5709,8 @@ Interval (in months) between two waves: for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ k=k+1; - /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/ + /* to clean */ + printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,codtab[cptcod][cptcov],nbcode); fprintf(ficrespl,"\n#******"); printf("\n#******"); fprintf(ficlog,"\n#******"); @@ -5979,6 +6022,7 @@ Interval (in months) between two waves: free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); } /* mle==-3 arrives here for freeing */ + endfree: free_matrix(prlim,1,nlstate,1,nlstate); free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); @@ -5996,7 +6040,6 @@ Interval (in months) between two waves: free_ivector(Tprod,1,15); free_ivector(Tvaraff,1,15); free_ivector(Tage,1,15); - free_ivector(Tcode,1,100); free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX); free_imatrix(codtab,1,100,1,10);