From 4d1da8bceb80178d7d5ca12ad6b8f3c1a6471c7f Mon Sep 17 00:00:00 2001 From: "N. Brouard" Date: Wed, 25 Jan 2006 00:51:50 +0000 Subject: [PATCH] (Module): Lots of cleaning and bugs added (Gompertz) --- src/imach.c | 212 +++++++++++++++++++++++++++------------------------- 1 file changed, 109 insertions(+), 103 deletions(-) diff --git a/src/imach.c b/src/imach.c index df7c6ff..1a5a087 100644 --- a/src/imach.c +++ b/src/imach.c @@ -1,6 +1,9 @@ /* $Id$ $State$ $Log$ + Revision 1.109 2006/01/24 19:37:15 brouard + (Module): Comments (lines starting with a #) are allowed in data. + Revision 1.108 2006/01/19 18:05:42 lievre Gnuplot problem appeared... To be fixed @@ -295,6 +298,7 @@ int popbased=0; int *wav; /* Number of waves for this individuual 0 is possible */ int maxwav; /* Maxim number of waves */ int jmin, jmax; /* min, max spacing between 2 waves */ +int ijmin, ijmax; /* Individuals having jmin and jmax */ int gipmx, gsw; /* Global variables on the number of contributions to the likelihood and the sum of weights (done by funcone)*/ int mle, weightopt; @@ -2225,8 +2229,14 @@ void concatwav(int wav[], int **dh, int **bh, int **mw, int **s, double *agedc fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); } k=k+1; - if (j >= jmax) jmax=j; - if (j <= jmin) jmin=j; + if (j >= jmax){ + jmax=j; + ijmax=i; + } + if (j <= jmin){ + jmin=j; + ijmin=i; + } sum=sum+j; /*if (j<0) printf("j=%d num=%d \n",j,i);*/ /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ @@ -2237,8 +2247,14 @@ void concatwav(int wav[], int **dh, int **bh, int **mw, int **s, double *agedc /* if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */ k=k+1; - if (j >= jmax) jmax=j; - else if (j <= jmin)jmin=j; + if (j >= jmax) { + jmax=j; + ijmax=i; + } + else if (j <= jmin){ + jmin=j; + ijmin=i; + } /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/ if(j<0){ @@ -2281,8 +2297,8 @@ void concatwav(int wav[], int **dh, int **bh, int **mw, int **s, double *agedc } /* end wave */ } jmean=sum/k; - printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean); - fprintf(ficlog,"Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean); + printf("Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, num[ijmin], jmax, num[ijmax], jmean); + fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean); } /*********** Tricode ****************************/ @@ -3993,6 +4009,7 @@ double gompertz(double x[]) { double A,B,L=0.0,sump=0.,num=0.; 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;*/ @@ -4005,14 +4022,15 @@ double gompertz(double x[]) for (i=1;i<=imx ; i++) { - if (cens[i]==1 & wav[i]>1) + if (cens[i] == 1 && wav[i]>1) A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); - if (cens[i]==0 & wav[i]>1) + if (cens[i] == 0 && wav[i]>1) A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))) +log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM); - if (wav[i]>1 & agecens[i]>15) { + /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */ + if (wav[i] > 1 ) { /* ??? */ L=L+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]);*/ } @@ -4081,7 +4099,7 @@ int main(int argc, char *argv[]) { 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; + int linei, month, year,iout; int jj, ll, li, lj, lk, imk; int numlinepar=0; /* Current linenumber of parameter file */ int itimes; @@ -4489,14 +4507,16 @@ int main(int argc, char *argv[]) i=1; linei=0; - while ((fgets(line, MAXLINE, fic) != NULL) ||((i >= firstobs) && (i <=lastobs))) { + while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) { linei=linei+1; - printf("IIIII= %d linei=%d\n",i,linei); 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);j--){;};line[j+1]=0; /* Trims blanks at end of line */ + 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); @@ -4513,66 +4533,48 @@ int main(int argc, char *argv[]) } s[j][i]=lval; - strcpy(line,stra); - 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\nShould be a year of exam at wave %d. Exiting.\n",lval, i,line,linei,j); - exit(1); - } - anint[j][i]=(double)(lval); - strcpy(line,stra); 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\nShould be a month of exam at wave %d. Exiting.\n",lval, i,line, linei,j); + if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ + } + else if(iout=sscanf(strb,".") != 0){ + month=99; + year=9999; + }else{ + printf("Error reading data around '%s'.at line number %ld %s for individual %d\nShould be a year of exam at wave %d. Exiting.\n",strb, i,line,linei,j); exit(1); } - mint[j][i]=(double)(lval); + anint[j][i]= (double) year; + mint[j][i]= (double)month; strcpy(line,stra); - } + } /* ENd Waves */ - 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\nShould be a year of death. Exiting.\n",lval, i,line,linei); + cutv(stra, strb,line,' '); + if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ + } + else if(iout=sscanf(strb,".") != 0){ + month=99; + year=9999; + }else{ + printf("Error reading data around '%s'.at line number %ld %s for individual %d\nShould be a year of exam at wave %d. Exiting.\n",strb, i,line,linei,j); exit(1); } - andc[i]=(double)(lval); + andc[i]=(double) year; + moisdc[i]=(double) month; strcpy(line,stra); 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\nShould be a month of death. Exiting.\n",lval,i,line, linei); - exit(1); - } - moisdc[i]=(double)(lval); - - strcpy(line,stra); - 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\nShould be a year of birth. Exiting.\n",lval, i,line, linei); - exit(1); + if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ } - annais[i]=(double)(lval); - - strcpy(line,stra); - 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\nShould be a month of birth. Exiting.\n",lval,i,line,linei); + else if(iout=sscanf(strb,".") != 0){ + month=99; + year=9999; + }else{ + printf("Error reading data around '%s'.at line number %ld %s for individual %d\nShould be a year of exam at wave %d. Exiting.\n",strb, i,line,linei,j); exit(1); } - moisnais[i]=(double)(lval); + annais[i]=(double)(year); + moisnais[i]=(double)(month); strcpy(line,stra); cutv(stra, strb,line,' '); @@ -4593,7 +4595,7 @@ int main(int argc, char *argv[]) printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a covar (meaning 0 for the reference or 1). Exiting.\n",lval, i,line,linei); exit(1); } - if(lval <0 || lval >1){ + if(lval <-1 || lval >1){ printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a value of the %d covar (meaning 0 for the reference or 1. IMaCh does not build design variables, do it your self). Exiting.\n",lval,i,line,linei,j); exit(1); } @@ -4608,7 +4610,6 @@ int main(int argc, char *argv[]) } else num[i]=atol(stra); - printf ("num [i] %ld %d\n",i, num[i]);fflush(stdout); /*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;}*/ @@ -4937,6 +4938,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ p=param[1][1]; /* *(*(*(param +1)+1)+0) */ globpr=0; /* To get the number ipmx of contributions and the sum of weights*/ + if (mle==-3){ ximort=matrix(1,NDIM,1,NDIM); cens=ivector(1,n); @@ -4946,9 +4948,9 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ for (i=1; i<=imx; i++){ dcwave[i]=-1; - for (j=1; j<=lastpass; j++) - if (s[j][i]>nlstate) { - dcwave[i]=j; + for (m=firstpass; m<=lastpass; m++) + if (s[m][i]>nlstate) { + dcwave[i]=m; /* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/ break; } @@ -4957,12 +4959,16 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ for (i=1; i<=imx; i++) { if (wav[i]>0){ ageexmed[i]=agev[mw[1][i]][i]; - j=wav[i];agecens[i]=1.; - if (ageexmed[i]>1 & wav[i]>0) agecens[i]=agev[mw[j][i]][i]; - cens[i]=1; - - if (ageexmed[i]<1) cens[i]=-1; - if (agedc[i]< AGESUP & agedc[i]>1 & dcwave[i]>firstpass & dcwave[i]<=lastpass) cens[i]=0 ; + j=wav[i]; + agecens[i]=1.; + + if (ageexmed[i]> 1 && wav[i] > 0){ + agecens[i]=agev[mw[j][i]][i]; + cens[i]= 1; + }else if (ageexmed[i]< 1) + cens[i]= -1; + if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass) + cens[i]=0 ; } else cens[i]=-1; } @@ -4971,29 +4977,29 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ for (j=1;j<=NDIM;j++) ximort[i][j]=(i == j ? 1.0 : 0.0); } - - p[1]=0.1; p[2]=0.1; + + p[1]=0.1; p[NDIM]=0.1; /*printf("%lf %lf", p[1], p[2]);*/ - printf("Powell\n"); fprintf(ficlog,"Powell\n"); - strcpy(filerespow,"pow-mort"); - strcat(filerespow,fileres); - if((ficrespow=fopen(filerespow,"w"))==NULL) { - printf("Problem with resultfile: %s\n", filerespow); - fprintf(ficlog,"Problem with resultfile: %s\n", filerespow); - } - fprintf(ficrespow,"# Powell\n# iter -2*LL"); - /* for (i=1;i<=nlstate;i++) - for(j=1;j<=nlstate+ndeath;j++) - if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); - */ - fprintf(ficrespow,"\n"); - + printf("Powell\n"); fprintf(ficlog,"Powell\n"); + strcpy(filerespow,"pow-mort"); + strcat(filerespow,fileres); + if((ficrespow=fopen(filerespow,"w"))==NULL) { + printf("Problem with resultfile: %s\n", filerespow); + fprintf(ficlog,"Problem with resultfile: %s\n", filerespow); + } + fprintf(ficrespow,"# Powell\n# iter -2*LL"); + /* for (i=1;i<=nlstate;i++) + for(j=1;j<=nlstate+ndeath;j++) + if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); + */ + fprintf(ficrespow,"\n"); + powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); fclose(ficrespow); - hesscov(matcov, p, NDIM,delti, 1e-4, gompertz); + hesscov(matcov, p, NDIM, delti, 1e-4, gompertz); for(i=1; i <=NDIM; i++) for(j=i+1;j<=NDIM;j++) @@ -5011,48 +5017,48 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ for (i=1;i<=NDIM;i++) printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); -lsurv=vector(1,AGESUP); + lsurv=vector(1,AGESUP); lpop=vector(1,AGESUP); tpop=vector(1,AGESUP); lsurv[agegomp]=100000; - - for (k=agegomp;k<=AGESUP;k++) { + + for (k=agegomp;k<=AGESUP;k++) { agemortsup=k; if (p[1]*exp(p[2]*(k-agegomp))>1) break; } - - for (k=agegomp;k=1 */ likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ -- 2.43.0