--- imach/src/imach.c 2002/02/26 17:11:54 1.25 +++ imach/src/imach.c 2002/02/28 17:49:07 1.27 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.25 2002/02/26 17:11:54 lievre Exp $ +/* $Id: imach.c,v 1.27 2002/02/28 17:49:07 lievre Exp $ Interpolated Markov Chain Short summary of the programme: @@ -94,7 +94,7 @@ int **dh; /* dh[mi][i] is number of step double jmean; /* Mean space between 2 waves */ double **oldm, **newm, **savm; /* Working pointers to matrices */ double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ -FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf; +FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; FILE *ficgp,*ficresprob,*ficpop; FILE *ficreseij; char filerese[FILENAMELENGTH]; @@ -1176,7 +1176,7 @@ void lubksb(double **a, int n, int *indx } /************ Frequencies ********************/ -void freqsummary(char fileres[], 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) +void freqsummary(char fileres[], 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 jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2) { /* Some frequencies */ int i, m, jk, k1,i1, j1, bool, z1,z2,j; @@ -1236,6 +1236,9 @@ void freqsummary(char fileres[], int ag } } } + + fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); + if (cptcovn>0) { fprintf(ficresp, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); @@ -1620,7 +1623,7 @@ void varevsij(char fileres[], double *** 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++) @@ -1632,7 +1635,7 @@ 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]; @@ -1684,7 +1687,7 @@ void varevsij(char fileres[], double *** free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar); free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); } /* End age */ - + free_vector(xp,1,npar); free_matrix(doldm,1,nlstate,1,npar); free_matrix(dnewm,1,nlstate,1,nlstate); @@ -1876,7 +1879,7 @@ if (i== 4) fprintf(ficresprob,"%.3e %.3e } free_vector(xp,1,npar); fclose(ficresprob); - exit(0); + } /******************* Printing html file ***********/ @@ -2143,6 +2146,261 @@ void movingaverage(double agedeb, double } + +/************** Forecasting ******************/ +prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double agemin, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){ + + int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h; + int *popage; + double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; + double *popeffectif,*popcount; + double ***p3mat; + char fileresf[FILENAMELENGTH]; + + agelim=AGESUP; +calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; + + prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); + + + strcpy(fileresf,"f"); + strcat(fileresf,fileres); + if((ficresf=fopen(fileresf,"w"))==NULL) { + printf("Problem with forecast resultfile: %s\n", fileresf); + } + printf("Computing forecasting: result on file '%s' \n", fileresf); + + if (cptcoveff==0) ncodemax[cptcoveff]=1; + + if (mobilav==1) { + mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + movingaverage(agedeb, fage, agemin, mobaverage); + } + + stepsize=(int) (stepm+YEARM-1)/YEARM; + if (stepm<=12) stepsize=1; + + agelim=AGESUP; + + hstepm=1; + hstepm=hstepm/stepm; + yp1=modf(dateintmean,&yp); + anprojmean=yp; + yp2=modf((yp1*12),&yp); + mprojmean=yp; + yp1=modf((yp2*30.5),&yp); + jprojmean=yp; + if(jprojmean==0) jprojmean=1; + if(mprojmean==0) jprojmean=1; + + fprintf(ficresf,"# Estimated date of observed prevalence: %.lf/%.lf/%.lf ",jprojmean,mprojmean,anprojmean); + + for(cptcov=1;cptcov<=i2;cptcov++){ + for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ + k=k+1; + fprintf(ficresf,"\n#******"); + for(j=1;j<=cptcoveff;j++) { + fprintf(ficresf," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + } + fprintf(ficresf,"******\n"); + fprintf(ficresf,"# StartingAge FinalAge"); + for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j); + + + for (cpt=0; cpt<=(anproj2-anproj1);cpt++) { + fprintf(ficresf,"\n"); + fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt); + + for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ + nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); + nhstepm = nhstepm/hstepm; + + p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + oldm=oldms;savm=savms; + hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); + + for (h=0; h<=nhstepm; h++){ + if (h==(int) (calagedate+YEARM*cpt)) { + fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm); + } + for(j=1; j<=nlstate+ndeath;j++) { + kk1=0.;kk2=0; + for(i=1; i<=nlstate;i++) { + if (mobilav==1) + kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod]; + else { + kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod]; + } + + } + if (h==(int)(calagedate+12*cpt)){ + fprintf(ficresf," %.3f", kk1); + + } + } + } + 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); +} +/************** Forecasting ******************/ +populforecast(char fileres[], double anproj1,double mproj1,double jproj1,double agemin, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){ + + int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h; + int *popage; + double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; + double *popeffectif,*popcount; + double ***p3mat,***tabpop,***tabpopprev; + char filerespop[FILENAMELENGTH]; + +tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); +tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + agelim=AGESUP; +calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; + + prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); + + + strcpy(filerespop,"pop"); + strcat(filerespop,fileres); + if((ficrespop=fopen(filerespop,"w"))==NULL) { + printf("Problem with forecast resultfile: %s\n", filerespop); + } + printf("Computing forecasting: result on file '%s' \n", filerespop); + + if (cptcoveff==0) ncodemax[cptcoveff]=1; + + if (mobilav==1) { + mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + movingaverage(agedeb, fage, agemin, mobaverage); + } + + stepsize=(int) (stepm+YEARM-1)/YEARM; + if (stepm<=12) stepsize=1; + + agelim=AGESUP; + + hstepm=1; + hstepm=hstepm/stepm; + + if (popforecast==1) { + if((ficpop=fopen(popfile,"r"))==NULL) { + printf("Problem with population file : %s\n",popfile);exit(0); + } + popage=ivector(0,AGESUP); + popeffectif=vector(0,AGESUP); + popcount=vector(0,AGESUP); + + i=1; + while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF) i=i+1; + + imx=i; + for (i=1; i=(agemin-((int)calagedate %12)/12.); agedeb--){ + nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); + nhstepm = nhstepm/hstepm; + + p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + oldm=oldms;savm=savms; + hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); + + for (h=0; h<=nhstepm; h++){ + if (h==(int) (calagedate+YEARM*cpt)) { + fprintf(ficrespop,"\n %.f ",agedeb+h*hstepm/YEARM*stepm); + } + for(j=1; j<=nlstate+ndeath;j++) { + kk1=0.;kk2=0; + for(i=1; i<=nlstate;i++) { + if (mobilav==1) + kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod]; + else { + kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod]; + } + } + if (h==(int)(calagedate+12*cpt)){ + tabpop[(int)(agedeb)][j][cptcod]=kk1; + /*fprintf(ficrespop," %.3f", kk1); + if (popforecast==1) fprintf(ficrespop," [%.f]", kk1*popeffectif[(int)agedeb+1]);*/ + } + } + for(i=1; i<=nlstate;i++){ + kk1=0.; + for(j=1; j<=nlstate;j++){ + kk1= kk1+Tabpop[(int)(agedeb)][j][cptcod]; + } + tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedate+12*cpt)*hstepm/YEARM*stepm-1)]; + } + +if (h==(int)(calagedate+12*cpt)) for(j=1; j<=nlstate;j++) fprintf(ficrespop," %.3f",tabpopprev[(int)(agedeb+1)][j][cptcod]); + } + free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + } + } + + /******/ + + for (cpt=1; cpt<=4;cpt++) { + fprintf(ficrespop,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt); + for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ + nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); + nhstepm = nhstepm/hstepm; + + p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + oldm=oldms;savm=savms; + hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); + for (h=0; h<=nhstepm; h++){ + if (h==(int) (calagedate+YEARM*cpt)) { + fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm); + } + for(j=1; j<=nlstate+ndeath;j++) { + kk1=0.;kk2=0; + for(i=1; i<=nlstate;i++) { + kk1=kk1+p3mat[i][j][h]*tabpopprev[(int)agedeb+1][i][cptcod]; + } + if (h==(int)(calagedate+12*cpt)) fprintf(ficresf," %.3f", kk1); + } + } + free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + } + } + } + } + + if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + + if (popforecast==1) { + free_ivector(popage,0,AGESUP); + free_vector(popeffectif,0,AGESUP); + free_vector(popcount,0,AGESUP); + } + free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + fclose(ficrespop); +} + /***********************************************/ /**************** Main Program *****************/ /***********************************************/ @@ -2165,7 +2423,7 @@ int main(int argc, char *argv[]) char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH]; char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH]; - char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH]; + char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH]; char filerest[FILENAMELENGTH]; char fileregp[FILENAMELENGTH]; @@ -2179,7 +2437,6 @@ int main(int argc, char *argv[]) int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; int mobilav=0,popforecast=0; int hstepm, nhstepm; - int *popage;/*boolprev=0 if date and zero if wave*/ double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2; double bage, fage, age, agelim, agebase; @@ -2195,9 +2452,8 @@ int main(int argc, char *argv[]) double **varpl; /* Variances of prevalence limits by age */ double *epj, vepp; double kk1, kk2; - double *popeffectif,*popcount; - double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,jprojmean,mprojmean,anprojmean, calagedate; - double yp,yp1,yp2; + double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2; + char version[80]="Imach version 0.7, February 2002, INED-EUROREVES "; char *alph[]={"a","a","b","c","d","e"}, str[4]; @@ -2763,11 +3019,11 @@ printf("Total number of individuals= %d, fputs(line,ficparo); } ungetc(c,ficpar); - fscanf(ficpar,"popforecast=%d popfile=%s starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mob_average=%d\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilav); -fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mob_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav); -fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mob_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav); + fscanf(ficpar,"popforecast=%d popfile=%s starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mov_average=%d\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilav); +fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav); +fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav); - freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2); + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); /*------------ gnuplot -------------*/ printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, agemin,agemax,fage, pathc,p); @@ -2899,142 +3155,20 @@ fprintf(ficres,"popforecast=%d popfile=% fclose(ficrespij); - if(stepm == 1) { - /*---------- Forecasting ------------------*/ - calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; - - /*printf("calage= %f", calagedate);*/ - - prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); - - - strcpy(fileresf,"f"); - strcat(fileresf,fileres); - if((ficresf=fopen(fileresf,"w"))==NULL) { - printf("Problem with forecast resultfile: %s\n", fileresf);goto end; - } - printf("Computing forecasting: result on file '%s' \n", fileresf); - - free_matrix(mint,1,maxwav,1,n); - free_matrix(anint,1,maxwav,1,n); - free_matrix(agev,1,maxwav,1,imx); - - /* Mobile average */ - - if (cptcoveff==0) ncodemax[cptcoveff]=1; - - if (mobilav==1) { - mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - movingaverage(agedeb, fage, agemin, mobaverage); - } - - stepsize=(int) (stepm+YEARM-1)/YEARM; - if (stepm<=12) stepsize=1; - - agelim=AGESUP; - /*hstepm=stepsize*YEARM; *//* Every year of age */ - hstepm=1; - hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */ - yp1=modf(dateintmean,&yp); - anprojmean=yp; - yp2=modf((yp1*12),&yp); - mprojmean=yp; - yp1=modf((yp2*30.5),&yp); - jprojmean=yp; - if(jprojmean==0) jprojmean=1; - if(mprojmean==0) jprojmean=1; - - fprintf(ficresf,"# Estimated date of observed prevalence: %.lf/%.lf/%.lf ",jprojmean,mprojmean,anprojmean); - - if (popforecast==1) { - if((ficpop=fopen(popfile,"r"))==NULL) { - printf("Problem with population file : %s\n",popfile);goto end; - } - popage=ivector(0,AGESUP); - popeffectif=vector(0,AGESUP); - popcount=vector(0,AGESUP); - - i=1; - while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF) - { - i=i+1; - } - imx=i; - - for (i=1; i=(agemin-((int)calagedate %12)/12.); agedeb--){ - nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); - nhstepm = nhstepm/hstepm; - - p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - oldm=oldms;savm=savms; - hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); - - for (h=0; h<=nhstepm; h++){ - if (h==(int) (calagedate+YEARM*cpt)) { - fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm); - } - for(j=1; j<=nlstate+ndeath;j++) { - kk1=0.;kk2=0; - for(i=1; i<=nlstate;i++) { - if (mobilav==1) - kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod]; - else { - kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod]; - - } - - if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb]; - } - - if (h==(int)(calagedate+12*cpt)){ - fprintf(ficresf," %.3f", kk1); - - if (popforecast==1) fprintf(ficresf," [%.f]", kk2); - } - } - } - free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - } - } - } - } - if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - if (popforecast==1) { - free_ivector(popage,0,AGESUP); - free_vector(popeffectif,0,AGESUP); - free_vector(popcount,0,AGESUP); - } - free_imatrix(s,1,maxwav+1,1,n); - free_vector(weight,1,n); - fclose(ficresf); - }/* End forecasting */ + /*---------- Forecasting ------------------*/ + if(stepm == 1) { + prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1); +populforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1); + free_matrix(mint,1,maxwav,1,n); + free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n); + free_vector(weight,1,n);} else{ erreur=108; printf("Error %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d\n", erreur, stepm); - } + - /*---------- Health expectancies and variances ------------*/ strcpy(filerest,"t"); @@ -3083,12 +3217,14 @@ fprintf(ficres,"popforecast=%d popfile=% evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k); vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); oldm=oldms;savm=savms; - varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k); - + varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k); + + + fprintf(ficrest,"#Total LEs with variances: e.. (std) "); for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); fprintf(ficrest,"\n"); - + hf=1; if (stepm >= YEARM) hf=stepm/YEARM; epj=vector(1,nlstate+1); @@ -3117,20 +3253,16 @@ fprintf(ficres,"popforecast=%d popfile=% } } } - - - - fclose(ficreseij); - fclose(ficresvij); + fclose(ficreseij); + fclose(ficresvij); fclose(ficrest); fclose(ficpar); free_vector(epj,1,nlstate+1); - /* scanf("%d ",i); */ - + /*------- Variance limit prevalence------*/ -strcpy(fileresvpl,"vpl"); + strcpy(fileresvpl,"vpl"); strcat(fileresvpl,fileres); if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { printf("Problem with variance prev lim resultfile: %s\n", fileresvpl); @@ -3138,19 +3270,19 @@ strcpy(fileresvpl,"vpl"); } printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl); - k=0; - for(cptcov=1;cptcov<=i1;cptcov++){ - for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ - k=k+1; - fprintf(ficresvpl,"\n#****** "); - for(j=1;j<=cptcoveff;j++) - fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); - fprintf(ficresvpl,"******\n"); - - varpl=matrix(1,nlstate,(int) bage, (int) fage); - oldm=oldms;savm=savms; + k=0; + for(cptcov=1;cptcov<=i1;cptcov++){ + for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ + k=k+1; + fprintf(ficresvpl,"\n#****** "); + for(j=1;j<=cptcoveff;j++) + fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficresvpl,"******\n"); + + varpl=matrix(1,nlstate,(int) bage, (int) fage); + oldm=oldms;savm=savms; varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k); - } + } } fclose(ficresvpl); @@ -3169,7 +3301,7 @@ strcpy(fileresvpl,"vpl"); free_matrix(matcov,1,npar,1,npar); free_vector(delti,1,npar); - + free_matrix(agev,1,maxwav,1,imx); free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); if(erreur >0)