--- imach/src/imach.c 2002/02/28 17:49:07 1.27 +++ imach/src/imach.c 2002/03/13 17:19:16 1.34 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.27 2002/02/28 17:49:07 lievre Exp $ +/* $Id: imach.c,v 1.34 2002/03/13 17:19:16 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -82,7 +82,7 @@ int cptcovn, cptcovage=0, cptcoveff=0,cp int npar=NPARMAX; int nlstate=2; /* Number of live states */ int ndeath=1; /* Number of dead states */ -int ncovmodel, ncov; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */ +int ncovmodel, ncovcol; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */ int popbased=0; int *wav; /* Number of waves for this individuual 0 is possible */ @@ -1185,7 +1185,7 @@ void freqsummary(char fileres[], int ag double pos, k2, dateintsum=0,k2cpt=0; FILE *ficresp; char fileresp[FILENAMELENGTH]; - + pp=vector(1,nlstate); probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); strcpy(fileresp,"p"); @@ -1307,7 +1307,7 @@ void freqsummary(char fileres[], int ag } /************ Prevalence ********************/ -void prevalence(int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, double calagedate) +void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, double calagedate) { /* Some frequencies */ int i, m, jk, k1, i1, j1, bool, z1,z2,j; @@ -1332,14 +1332,14 @@ void prevalence(int agemin, int agemax, for (jk=-1; jk<=nlstate+ndeath; jk++) for(m=agemin; m <= agemax+3; m++) freq[i][jk][m]=0; - + for (i=1; i<=imx; i++) { bool=1; if (cptcovn>0) { for (z1=1; z1<=cptcoveff; z1++) if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) bool=0; - } + } if (bool==1) { for(m=firstpass; m<=lastpass; m++){ k2=anint[m][i]+(mint[m][i]/12.); @@ -1347,12 +1347,11 @@ void prevalence(int agemin, int agemax, if(agev[m][i]==0) agev[m][i]=agemax+1; if(agev[m][i]==1) agev[m][i]=agemax+2; freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i]; - freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i]; + /* freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i]; */ } } } } - for(i=(int)agemin; i <= (int)agemax+3; i++){ for(jk=1; jk <=nlstate ; jk++){ for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) @@ -1444,7 +1443,7 @@ void concatwav(int wav[], int **dh, int if (j >= jmax) jmax=j; if (j <= jmin) jmin=j; sum=sum+j; - /* if (j<10) printf("j=%d num=%d ",j,i); */ + /*if (j<0) printf("j=%d num=%d \n",j,i); */ } } else{ @@ -1452,7 +1451,7 @@ void concatwav(int wav[], int **dh, int k=k+1; if (j >= jmax) jmax=j; else if (j <= jmin)jmin=j; - /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ + /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ sum=sum+j; } jk= j/stepm; @@ -1514,7 +1513,7 @@ void tricode(int *Tvar, int **nbcode, in ij=1; for (i=1; i<=10; i++) { - if((Ndum[i]!=0) && (i<=ncov)){ + if((Ndum[i]!=0) && (i<=ncovcol)){ Tvaraff[ij]=i; ij++; } @@ -1528,7 +1527,7 @@ void tricode(int *Tvar, int **nbcode, in void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij) { /* Health expectancies */ - int i, j, nhstepm, hstepm, h; + int i, j, nhstepm, hstepm, h, nstepm, k; double age, agelim,hf; double ***p3mat; @@ -1539,34 +1538,44 @@ void evsij(char fileres[], double ***eij fprintf(ficreseij," %1d-%1d",i,j); fprintf(ficreseij,"\n"); - hstepm=1*YEARM; /* Every j years of age (in month) */ - hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */ + k=1; /* For example stepm=6 months */ + hstepm=k*YEARM; /* (a) Every k years of age (in months), for example every k=2 years 24 m */ + hstepm=stepm; /* or (b) We decided to compute the life expectancy with the smallest unit */ + /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. + nhstepm is the number of hstepm from age to agelim + nstepm is the number of stepm from age to agelin. + Look at hpijx to understand the reason of that which relies in memory size + and note for a fixed period like k years */ + /* We decided (b) to get a life expectancy respecting the most precise curvature of the + survival function given by stepm (the optimization length). Unfortunately it + means that if the survival funtion is printed only each two years of age and if + you sum them up and add 1 year (area under the trapezoids) you won't get the same + results. So we changed our mind and took the option of the best precision. + */ + hstepm=hstepm/stepm; /* Typically in stepm units, if k= 2 years, = 2/6 months = 4 */ agelim=AGESUP; for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ /* nhstepm age range expressed in number of stepm */ - nhstepm=(int) rint((agelim-age)*YEARM/stepm); - /* Typically if 20 years = 20*12/6=40 stepm */ - if (stepm >= YEARM) hstepm=1; - nhstepm = nhstepm/hstepm;/* Expressed in hstepm, typically 40/4=10 */ + nstepm=(int) rint((agelim-age)*YEARM/stepm); + /* Typically if 20 years nstepm = 20*12/6=40 stepm */ + /* if (stepm >= YEARM) hstepm=1;*/ + nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); /* Computed by stepm unit matrices, product of hstepm matrices, stored in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij); - - + hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) - for (h=0, eij[i][j][(int)age]=0; h<=nhstepm; h++){ - eij[i][j][(int)age] +=p3mat[i][j][h]; + for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ + eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf; + /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ } - - hf=1; - if (stepm >= YEARM) hf=stepm/YEARM; - fprintf(ficreseij,"%.0f",age ); + fprintf(ficreseij,"%3.0f",age ); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++){ - fprintf(ficreseij," %.4f", hf*eij[i][j][(int)age]); + fprintf(ficreseij," %9.4f", eij[i][j][(int)age]); } fprintf(ficreseij,"\n"); free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); @@ -1883,7 +1892,7 @@ fclose(ficresprob); } /******************* Printing html file ***********/ -void printinghtml(char fileres[], char title[], char datafile[], int firstpass, int lastpass, int stepm, int weightopt, char model[],int imx,int jmin, int jmax, double jmeanint,char optionfile[],char optionfilehtm[] ){ +void printinghtml(char fileres[], char title[], char datafile[], int firstpass, int lastpass, int stepm, int weightopt, char model[],int imx,int jmin, int jmax, double jmeanint,char optionfile[],char optionfilehtm[],char rfileres[] ){ int jj1, k1, i1, cpt; FILE *fichtm; /*char optionfilehtm[FILENAMELENGTH];*/ @@ -1894,7 +1903,7 @@ void printinghtml(char fileres[], char t printf("Problem with %s \n",optionfilehtm), exit(0); } - fprintf(fichtm,"