--- imach/src/imach.c 2003/03/28 13:33:56 1.72 +++ imach/src/imach.c 2003/04/08 14:06:50 1.73 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.72 2003/03/28 13:33:56 brouard Exp $ +/* $Id: imach.c,v 1.73 2003/04/08 14:06:50 lievre Exp $ Interpolated Markov Chain Short summary of the programme: @@ -83,7 +83,7 @@ #define ODIRSEPARATOR '\\' #endif -char version[80]="Imach version 0.93, February 2003, INED-EUROREVES "; +char version[80]="Imach version 0.94, February 2003, INED-EUROREVES "; int erreur; /* Error number */ int nvar; int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; @@ -1429,12 +1429,13 @@ void freqsummary(char fileres[], int ag int i, m, jk, k1,i1, j1, bool, z1,z2,j; int first; double ***freq; /* Frequencies */ - double *pp; - double pos, k2, dateintsum=0,k2cpt=0; + double *pp, **prop; + double pos,posprop, k2, dateintsum=0,k2cpt=0; FILE *ficresp; char fileresp[FILENAMELENGTH]; pp=vector(1,nlstate); + prop=matrix(1,nlstate,agemin,agemax+3); probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); strcpy(fileresp,"p"); strcat(fileresp,fileres); @@ -1460,6 +1461,10 @@ void freqsummary(char fileres[], int ag for (jk=-1; jk<=nlstate+ndeath; jk++) for(m=agemin; m <= agemax+3; m++) freq[i][jk][m]=0; + + for (i=1; i<=nlstate; i++) + for(m=agemin; m <= agemax+3; m++) + prop[i][m]=0; dateintsum=0; k2cpt=0; @@ -1476,6 +1481,7 @@ void freqsummary(char fileres[], int ag if ((k2>=dateprev1) && (k2<=dateprev2)) { if(agev[m][i]==0) agev[m][i]=agemax+1; if(agev[m][i]==1) agev[m][i]=agemax+2; + if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i]; if (m=1.e-5){ if(first==1) @@ -1549,12 +1556,12 @@ void freqsummary(char fileres[], int ag } if( i <= (int) agemax){ if(pos>=1.e-5){ - fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos); + fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); probs[i][jk][j1]= pp[jk]/pos; /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/ } else - fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos); + fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop); } } @@ -1578,7 +1585,7 @@ void freqsummary(char fileres[], int ag fclose(ficresp); free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); free_vector(pp,1,nlstate); - + free_matrix(prop,1,nlstate,(int) agemin,(int) agemax+3); /* End of Freq */ } @@ -1592,12 +1599,12 @@ void prevalence(int agemin, float agemax int i, m, jk, k1, i1, j1, bool, z1,z2,j; double ***freq; /* Frequencies */ - double *pp; - double pos; + double *pp, **prop; + double pos,posprop; double y2; /* in fractional years */ pp=vector(1,nlstate); - + prop=matrix(1,nlstate,agemin,agemax+3); freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3); j1=0; @@ -1608,10 +1615,9 @@ void prevalence(int agemin, float agemax for(i1=1; i1<=ncodemax[k1];i1++){ j1++; - for (i=-1; i<=nlstate+ndeath; i++) - for (jk=-1; jk<=nlstate+ndeath; jk++) - for(m=agemin; m <= agemax+3; m++) - freq[i][jk][m]=0; + for (i=1; i<=nlstate; i++) + for(m=agemin; m <= agemax+3; m++) + prop[i][m]=0; for (i=1; i<=imx; i++) { /* Each individual */ bool=1; @@ -1626,35 +1632,24 @@ void prevalence(int agemin, float agemax if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ if(agev[m][i]==0) agev[m][i]=agemax+1; if(agev[m][i]==1) agev[m][i]=agemax+2; - if (m0 && s[m][i]<=nlstate) { + prop[s[m][i]][(int)agev[m][i]] += weight[i]; + prop[s[m][i]][(int)(agemax+3)] += weight[i]; } } } /* end selection of waves */ } } 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++) - pp[jk] += freq[jk][m][i]; - } - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pos=0; m <=0 ; m++) - pos += freq[jk][m][i]; - } - for(jk=1; jk <=nlstate ; jk++){ - for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++) - pp[jk] += freq[jk][m][i]; + for(jk=1,posprop=0; jk <=nlstate ; jk++) { + posprop += prop[jk][i]; } - for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk]; - for(jk=1; jk <=nlstate ; jk++){ if( i <= (int) agemax){ - if(pos>=1.e-5){ - probs[i][jk][j1]= pp[jk]/pos; + if(posprop>=1.e-5){ + probs[i][jk][j1]= prop[jk][i]/posprop; } } }/* end jk */ @@ -1665,7 +1660,7 @@ void prevalence(int agemin, float agemax free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); free_vector(pp,1,nlstate); - + free_matrix(prop,1,nlstate,(int) agemin,(int) agemax+3); } /* End of Freq */ /************* Waves Concatenation ***************/ @@ -1744,6 +1739,7 @@ void concatwav(int wav[], int **dh, int if (j >= jmax) jmax=j; else if (j <= jmin)jmin=j; /* 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]);*/ sum=sum+j; } jk= j/stepm; @@ -3048,7 +3044,7 @@ prevforecast(char fileres[], double anpr dateprev1 dateprev2 range of dates during which prevalence is computed anproj2 year of en of projection (same day and month as proj1). */ - int yearp, stepsize, hstepm, nhstepm, j, k, c, cptcod, i, h; + int yearp, stepsize, hstepm, nhstepm, j, k, c, cptcod, i, h, i1; int *popage; double agec; /* generic age */ double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; @@ -3093,11 +3089,15 @@ prevforecast(char fileres[], double anpr jprojmean=yp; if(jprojmean==0) jprojmean=1; if(mprojmean==0) jprojmean=1; + + i1=cptcoveff; + if (cptcovn < 1){i1=1;} fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); fprintf(ficresf,"#****** Routine prevforecast **\n"); - for(cptcov=1, k=0;cptcov<=cptcoveff;cptcov++){ + + for(cptcov=1, k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ k=k+1; fprintf(ficresf,"\n#******");