From 7d0c0715fddc718769ae1a14478f605efb229c50 Mon Sep 17 00:00:00 2001 From: "N. Brouard" Date: Sat, 20 Jun 2009 16:22:47 +0000 Subject: [PATCH] Some dimensions resccaled --- src/imach.c | 179 ++++++++++++++++++++++++++++++++-------------------- 1 file changed, 112 insertions(+), 67 deletions(-) diff --git a/src/imach.c b/src/imach.c index 737677b..07d193b 100644 --- a/src/imach.c +++ b/src/imach.c @@ -1,6 +1,11 @@ /* $Id$ $State$ $Log$ + 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 + V2+V3*age+V1+V4 strb=V3*age+V1+V4 working better. + Revision 1.129 2007/08/31 13:49:27 lievre Modification of the way of exiting when the covariate is not binary in order to see on the window the error message before exiting @@ -367,7 +372,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 @@ -392,7 +397,7 @@ extern int errno; /* $Id$ */ /* $State$ */ -char version[]="Imach version 0.98i, June 2006, INED-EUROREVES-Institut de longevite "; +char version[]="Imach version 0.98k, June 2006, INED-EUROREVES-Institut de longevite "; char fullversion[]="$Revision$ $Date$"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; @@ -516,7 +521,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 */ @@ -1166,7 +1171,7 @@ double **prevalim(double **prlim, int nlstate, double x[], double age, double ** 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 */ @@ -1248,10 +1253,14 @@ double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate ) 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; @@ -2470,9 +2480,9 @@ void tricode(int *Tvar, int **nbcode, int imx) 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]++; } @@ -2483,8 +2493,8 @@ void tricode(int *Tvar, int **nbcode, int imx) ij++; } } - - cptcoveff=ij-1; /*Number of simple covariates*/ + ij--; + cptcoveff=ij; /*Number of simple covariates*/ } /*********** Health Expectancies ****************/ @@ -3082,9 +3092,9 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double 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); */ @@ -3261,7 +3271,7 @@ void varprob(char optionfilefiname[], double **matcov, double x[], double delti[ 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); @@ -3414,7 +3424,7 @@ To be simple, these graphs help to understand the significativity of each parame /* 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); @@ -3676,7 +3686,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,subdirf2(fil 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 ++) { @@ -4626,7 +4636,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)*/ @@ -4752,6 +4768,9 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin 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) @@ -4815,7 +4834,7 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin 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); @@ -4839,12 +4858,17 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin 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; @@ -4856,8 +4880,9 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin 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; @@ -4871,8 +4896,9 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin 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; @@ -4886,7 +4912,8 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin 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); @@ -4897,18 +4924,25 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin 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 \ @@ -4920,6 +4954,15 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin 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); @@ -4958,7 +5001,7 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin 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); @@ -4974,7 +5017,7 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin 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; } @@ -5157,7 +5200,6 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin /* 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); @@ -5167,12 +5209,13 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin 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]]); } } @@ -5182,7 +5225,7 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin 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"); } @@ -5210,7 +5253,8 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin 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 */ @@ -5389,8 +5433,8 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ } /* 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]); @@ -5662,7 +5706,8 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ 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#******"); @@ -5974,6 +6019,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ 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); @@ -5991,7 +6037,6 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ 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); -- 2.43.0