--- imach096d/src/imach.c 2002/02/20 17:02:08 1.13 +++ imach096d/src/imach.c 2002/02/20 17:05:44 1.14 @@ -73,6 +73,7 @@ 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 popbased=0, fprev,lprev; int *wav; /* Number of waves for this individuual 0 is possible */ int maxwav; /* Maxim number of waves */ @@ -1191,7 +1192,7 @@ void freqsummary(char fileres[], int ag bool=0; } if (bool==1) { - for(m=firstpass; m<=lastpass-1; m++){ + for(m=fprev; m<=lprev; m++){ 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]] += weight[i]; @@ -1215,7 +1216,7 @@ void freqsummary(char fileres[], int ag printf("Age %d", i); for(jk=1; jk <=nlstate ; jk++){ for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) - pp[jk] += freq[jk][m][i]; + pp[jk] += freq[jk][m][i]; } for(jk=1; jk <=nlstate ; jk++){ for(m=-1, pos=0; m <=0 ; m++) @@ -1225,10 +1226,12 @@ void freqsummary(char fileres[], int ag else printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); } - for(jk=1; jk <=nlstate ; jk++){ - for(m=1, pp[jk]=0; m <=nlstate+ndeath; m++) + + 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,pos=0; jk <=nlstate ; jk++) pos += pp[jk]; for(jk=1; jk <=nlstate ; jk++){ @@ -1492,6 +1495,12 @@ void varevsij(char fileres[], double *** } hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); + + if (popbased==1) { + for(i=1; i<=nlstate;i++) + prlim[i][i]=probs[(int)age][i][ij]; + } + for(j=1; j<= nlstate; j++){ for(h=0; h<=nhstepm; h++){ for(i=1, gp[h][j]=0.;i<=nlstate;i++) @@ -1503,12 +1512,19 @@ void varevsij(char fileres[], double *** xp[i] = x[i] - (i==theta ?delti[theta]:0); hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); + + if (popbased==1) { + for(i=1; i<=nlstate;i++) + prlim[i][i]=probs[(int)age][i][ij]; + } + for(j=1; j<= nlstate; j++){ for(h=0; h<=nhstepm; h++){ for(i=1, gm[h][j]=0.;i<=nlstate;i++) gm[h][j] += prlim[i][i]*p3mat[i][j][h]; } } + for(j=1; j<= nlstate; j++) for(h=0; h<=nhstepm; h++){ gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; @@ -1773,9 +1789,10 @@ int main() int c, h , cpt,l; int ju,jl, mi; int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij; - int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; - + int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; + int mobilav=0, fprevfore=1, lprevfore=1; int hstepm, nhstepm; + double bage, fage, age, agelim, agebase; double ftolpl=FTOL; double **prlim; @@ -1852,7 +1869,25 @@ split(pathtot, path,optionfile); fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt,model); fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt,model); - +while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + puts(line); + fputs(line,ficparo); + } + ungetc(c,ficpar); + + fscanf(ficpar,"fprevalence=%d lprevalence=%d pop_based=%d\n",&fprev,&lprev,&popbased); + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + puts(line); + fputs(line,ficparo); + } + ungetc(c,ficpar); + + fscanf(ficpar,"fprevalence=%d lprevalence=%d mob_average=%d\n",&fprevfore,&lprevfore,&mobilav); + covar=matrix(0,NCOVMAX,1,n); cptcovn=0; if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; @@ -2014,8 +2049,8 @@ split(pathtot, path,optionfile); if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3; if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3; if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3; - } - for (i=1; i<=imx; i++) printf("%d %.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]));*/ + } + for (i=1; i<=imx; i++) printf("%d %.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]));*/ /* Calculation of the number of parameter from char model*/ Tvar=ivector(1,15); @@ -2525,7 +2560,9 @@ Interval (in months) between two waves: - Life expectancies by age and initial health status: e%s
- Variances of life expectancies by age and initial health status: v%s
- Health expectancies with their variances: t%s
- - Standard deviation of stationary prevalences: vpl%s

",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres); + - Standard deviation of stationary prevalences: vpl%s
+ - Prevalences forecasting: f%s
+
",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres); fprintf(fichtm,"
  • Graphs
  • "); @@ -2663,7 +2700,6 @@ fclose(fichtm); fclose(ficrespij); - exit(0); /*---------- Forecasting ------------------*/ strcpy(fileresf,"f"); @@ -2682,26 +2718,25 @@ fclose(fichtm); if (cptcoveff==0) ncodemax[cptcoveff]=1; - mobaverage= ma3x(1,130 ,1,8, 1,8); - for (agedeb=bage+3; agedeb<=fage-2; agedeb++) - for (i=1; i<=nlstate;i++) - for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++) - mobaverage[(int)agedeb][i][cptcod]=0.; - - for (agedeb=bage+4; agedeb<=fage; agedeb++){ - for (i=1; i<=nlstate;i++){ - for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ - for (cpt=0;cpt<=4;cpt++){ - mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod]; + if (mobilav==1) { + mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + for (agedeb=bage+3; agedeb<=fage-2; agedeb++) + for (i=1; i<=nlstate;i++) + for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++) + mobaverage[(int)agedeb][i][cptcod]=0.; + + for (agedeb=bage+4; agedeb<=fage; agedeb++){ + for (i=1; i<=nlstate;i++){ + for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ + for (cpt=0;cpt<=4;cpt++){ + mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod]; } mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5; + } } - } + } } -/* if (cptcod==2) printf("m=%f p=%f %d age=%d ",mobaverage[(int)agedeb-2][i][cptcod],probs[(int)agedeb-cpt][i][cptcod],cpt,(int)agedeb-2);*/ - - stepsize=(int) (stepm+YEARM-1)/YEARM; if (stepm<=24) stepsize=2; @@ -2725,12 +2760,18 @@ fclose(fichtm); for (agedeb=fage; agedeb>=bage; agedeb--){ fprintf(ficresf,"\n%d %.f %.f 0 ",k,agedeb, agedeb); + if (mobilav==1) { for(j=1; j<=nlstate;j++) - fprintf(ficresf,"%.3f ",mobaverage[(int)agedeb][j][cptcod]); - } + fprintf(ficresf,"%.5f ",mobaverage[(int)agedeb][j][cptcod]); + } + else { + for(j=1; j<=nlstate;j++) + fprintf(ficresf,"%.5f ",probs[(int)agedeb][j][cptcod]); + } + for(j=1; j<=ndeath;j++) fprintf(ficresf,"0."); - - for (cpt=1; cpt<=8;cpt++) + } + for (cpt=1; cpt<=NCOVMAX;cpt++) for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */ nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ @@ -2749,20 +2790,19 @@ fclose(fichtm); for(j=1; j<=nlstate+ndeath;j++) { kk1=0.; for(i=1; i<=nlstate;i++) { - /* kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];*/ - kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod]; - } - - if (h*hstepm/YEARM*stepm==cpt) - fprintf(ficresf," %.5f ", kk1); - } + if (mobilav==1) + kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod]; + else kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod]; + } + if (h*hstepm/YEARM*stepm==cpt) fprintf(ficresf," %.5f ", kk1); } - free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); } + free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); } } + } + if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); fclose(ficresf); - /*---------- Health expectancies and variances ------------*/ strcpy(filerest,"t"); @@ -2822,6 +2862,11 @@ fclose(fichtm); epj=vector(1,nlstate+1); for(age=bage; age <=fage ;age++){ prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); + if (popbased==1) { + for(i=1; i<=nlstate;i++) + prlim[i][i]=probs[(int)age][i][k]; + } + fprintf(ficrest," %.0f",age); for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ for(i=1, epj[j]=0.;i <=nlstate;i++) {