--- imach/src/imach.c 2002/07/23 23:59:37 1.53 +++ imach/src/imach.c 2002/07/26 12:29:55 1.58 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.53 2002/07/23 23:59:37 brouard Exp $ +/* $Id: imach.c,v 1.58 2002/07/26 12:29:55 lievre Exp $ Interpolated Markov Chain Short summary of the programme: @@ -39,7 +39,7 @@ hPijx. Also this programme outputs the covariance matrix of the parameters but also - of the life expectancies. It also computes the prevalence limits. + of the life expectancies. It also computes the stable prevalence. Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). Institut national d'études démographiques, Paris. @@ -60,7 +60,7 @@ /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ #define FILENAMELENGTH 80 /*#define DEBUG*/ -#define unix +#define windows #define GLOCK_ERROR_NOPATH -1 /* empty path */ #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ @@ -83,7 +83,7 @@ #define ODIRSEPARATOR '\\' #endif -char version[80]="Imach version 0.8j, July 2002, INED-EUROREVES "; +char version[80]="Imach version 0.8k, July 2002, INED-EUROREVES "; int erreur; /* Error number */ int nvar; int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; @@ -163,7 +163,7 @@ int estepm; int m,nb; int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage; double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint; -double **pmmij, ***probs, ***mobaverage; +double **pmmij, ***probs; double dateintmean=0; double *weight; @@ -701,12 +701,12 @@ void powell(double p[], double **xi, int printf("\n"); fprintf(ficlog,"\n"); #endif - } + } } } } -/**** Prevalence limit ****************/ +/**** Prevalence limit (stable prevalence) ****************/ double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij) { @@ -1284,7 +1284,7 @@ void freqsummary(char fileres[], int ag if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) bool=0; } - if (bool==1) { + if (bool==1){ for(m=firstpass; m<=lastpass; m++){ k2=anint[m][i]+(mint[m][i]/12.); if ((k2>=dateprev1) && (k2<=dateprev2)) { @@ -1575,31 +1575,35 @@ void concatwav(int wav[], int **dh, int /*********** Tricode ****************************/ void tricode(int *Tvar, int **nbcode, int imx) { - int Ndum[20],ij=1, k, j, i; + + int Ndum[20],ij=1, k, j, i, maxncov=19; int cptcode=0; cptcoveff=0; - for (k=0; k<19; k++) Ndum[k]=0; + for (k=0; k cptcode) cptcode=ij; + if (ij > cptcode) cptcode=ij; /* getting the maximum of covariable + Tvar[j]. If V=sex and male is 0 and + female is 1, then cptcode=1.*/ } for (i=0; i<=cptcode; i++) { - if(Ndum[i]!=0) ncodemax[j]++; + if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariates. In fact ncodemax[j]=2 (dichotom. variables) but it can be more */ } - ij=1; - + ij=1; for (i=1; i<=ncodemax[j]; i++) { - for (k=0; k<=19; k++) { + for (k=0; k<= maxncov; k++) { if (Ndum[k] != 0) { nbcode[Tvar[j]][ij]=k; + /* store the modality in an array. k is a modality. If we have model=V1+V1*sex then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ ij++; } @@ -1608,22 +1612,23 @@ void tricode(int *Tvar, int **nbcode, in } } - for (k=0; k<19; k++) Ndum[k]=0; + for (k=0; k< maxncov; k++) Ndum[k]=0; - for (i=1; i<=ncovmodel-2; i++) { + for (i=1; i<=ncovmodel-2; i++) { + /* Listing of all covariables in staement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ ij=Tvar[i]; - Ndum[ij]++; + Ndum[ij]++; } ij=1; - for (i=1; i<=10; i++) { + for (i=1; i<= maxncov; i++) { if((Ndum[i]!=0) && (i<=ncovcol)){ - Tvaraff[ij]=i; + Tvaraff[ij]=i; /*For printing */ ij++; } } - cptcoveff=ij-1; + cptcoveff=ij-1; /*Number of simple covariates*/ } /*********** Health Expectancies ****************/ @@ -1819,21 +1824,24 @@ void varevsij(char optionfilefiname[], d double ***mobaverage; int theta; char digit[4]; - char digitp[16]; + char digitp[25]; char fileresprobmorprev[FILENAMELENGTH]; - if(popbased==1) - strcpy(digitp,"-populbased-"); - else + if(popbased==1){ + if(mobilav!=0) + strcpy(digitp,"-populbased-mobilav-"); + else strcpy(digitp,"-populbased-nomobil-"); + } + else strcpy(digitp,"-stablbased-"); - if(mobilav==1) - strcat(digitp,"mobilav-"); - else - strcat(digitp,"nomobil-"); - if (mobilav==1) { + + if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - movingaverage(probs, bage, fage, mobaverage); + if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ + fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); + printf(" Error in movingaverage mobilav=%d\n",mobilav); + } } strcpy(fileresprobmorprev,"prmorprev"); @@ -1928,10 +1936,10 @@ void varevsij(char optionfilefiname[], d prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); if (popbased==1) { - if(mobilav !=1){ + if(mobilav ==0){ for(i=1; i<=nlstate;i++) prlim[i][i]=probs[(int)age][i][ij]; - }else{ /* mobilav=1 */ + }else{ /* mobilav */ for(i=1; i<=nlstate;i++) prlim[i][i]=mobaverage[(int)age][i][ij]; } @@ -1956,10 +1964,10 @@ void varevsij(char optionfilefiname[], d prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); if (popbased==1) { - if(mobilav !=1){ + if(mobilav ==0){ for(i=1; i<=nlstate;i++) prlim[i][i]=probs[(int)age][i][ij]; - }else{ /* mobilav=1 */ + }else{ /* mobilav */ for(i=1; i<=nlstate;i++) prlim[i][i]=mobaverage[(int)age][i][ij]; } @@ -2025,10 +2033,10 @@ void varevsij(char optionfilefiname[], d prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij); if (popbased==1) { - if(mobilav !=1){ + if(mobilav ==0){ for(i=1; i<=nlstate;i++) prlim[i][i]=probs[(int)age][i][ij]; - }else{ /* mobilav=1 */ + }else{ /* mobilav */ for(i=1; i<=nlstate;i++) prlim[i][i]=mobaverage[(int)age][i][ij]; } @@ -2084,11 +2092,10 @@ void varevsij(char optionfilefiname[], d free_matrix(doldmp,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); free_matrix(dnewmp,nlstate+1,nlstate+ndeath,1,npar); free_matrix(varppt,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); - free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); fclose(ficresprobmorprev); fclose(ficgp); fclose(fichtm); - } /************ Variance of prevlim ******************/ @@ -2106,7 +2113,7 @@ void varprevlim(char fileres[], double * double age,agelim; int theta; - fprintf(ficresvpl,"# Standard deviation of prevalence's limit\n"); + fprintf(ficresvpl,"# Standard deviation of stable prevalences \n"); fprintf(ficresvpl,"# Age"); for(i=1; i<=nlstate;i++) fprintf(ficresvpl," %1d-%1d",i,i); @@ -2599,9 +2606,9 @@ void printinggnuplot(char fileres[], dou fprintf(ficlog,"Problem with file %s",optionfilegnuplot); } -#ifdef windows + /*#ifdef windows */ fprintf(ficgp,"cd \"%s\" \n",pathc); -#endif + /*#endif */ m=pow(2,cptcoveff); /* 1eme*/ @@ -2614,7 +2621,7 @@ m=pow(2,cptcoveff); if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); else fprintf(ficgp," \%%*lf (\%%*lf)"); } - fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1); + fprintf(ficgp,"\" t\"Stable prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); else fprintf(ficgp," \%%*lf (\%%*lf)"); @@ -2768,27 +2775,45 @@ m=pow(2,cptcoveff); /*************** Moving average **************/ -void movingaverage(double ***probs, double bage,double fage, double ***mobaverage){ +int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav){ int i, cpt, cptcod; + int modcovmax =1; + int mobilavrange, mob; double age; - for (age=bage; age<=fage; age++) - for (i=1; i<=nlstate;i++) - for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++) - mobaverage[(int)age][i][cptcod]=0.; - - for (age=bage+4; age<=fage; age++){ - for (i=1; i<=nlstate;i++){ - for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ - for (cpt=0;cpt<=4;cpt++){ - mobaverage[(int)age-2][i][cptcod]=mobaverage[(int)age-2][i][cptcod]+probs[(int)age-cpt][i][cptcod]; + + modcovmax=2*cptcoveff;/* Max number of modalities. We suppose + a covariate has 2 modalities */ + if (cptcovn<1) modcovmax=1; /* At least 1 pass */ + + if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){ + if(mobilav==1) mobilavrange=5; /* default */ + else mobilavrange=mobilav; + for (age=bage; age<=fage; age++) + for (i=1; i<=nlstate;i++) + for (cptcod=1;cptcod<=modcovmax;cptcod++) + mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod]; + /* We keep the original values on the extreme ages bage, fage and for + fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2 + we use a 5 terms etc. until the borders are no more concerned. + */ + for (mob=3;mob <=mobilavrange;mob=mob+2){ + for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ + for (i=1; i<=nlstate;i++){ + for (cptcod=1;cptcod<=modcovmax;cptcod++){ + mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod]; + for (cpt=1;cpt<=(mob-1)/2;cpt++){ + mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod]; + mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod]; + } + mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob; + } } - mobaverage[(int)age-2][i][cptcod]=mobaverage[(int)age-2][i][cptcod]/5; - } - } - } - -} + }/* end age */ + }/* end mob */ + }else return -1; + return 0; +}/* End movingaverage */ /************** Forecasting ******************/ @@ -2799,10 +2824,11 @@ prevforecast(char fileres[], double anpr double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; double *popeffectif,*popcount; double ***p3mat; + double ***mobaverage; char fileresf[FILENAMELENGTH]; agelim=AGESUP; -calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; + calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); @@ -2818,9 +2844,12 @@ calagedate=(anproj1+mproj1/12.+jproj1/36 if (cptcoveff==0) ncodemax[cptcoveff]=1; - if (mobilav==1) { + if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - movingaverage(probs, ageminpar,fage, mobaverage); + if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ + fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); + printf(" Error in movingaverage mobilav=%d\n",mobilav); + } } stepsize=(int) (stepm+YEARM-1)/YEARM; @@ -2891,7 +2920,7 @@ calagedate=(anproj1+mproj1/12.+jproj1/36 } } - if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); fclose(ficresf); } @@ -2903,6 +2932,7 @@ populforecast(char fileres[], double anp double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; double *popeffectif,*popcount; double ***p3mat,***tabpop,***tabpopprev; + double ***mobaverage; char filerespop[FILENAMELENGTH]; tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); @@ -2924,9 +2954,12 @@ populforecast(char fileres[], double anp if (cptcoveff==0) ncodemax[cptcoveff]=1; - if (mobilav==1) { + if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - movingaverage(probs, ageminpar, fage, mobaverage); + if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ + fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); + printf(" Error in movingaverage mobilav=%d\n",mobilav); + } } stepsize=(int) (stepm+YEARM-1)/YEARM; @@ -3039,7 +3072,7 @@ populforecast(char fileres[], double anp } } - if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); if (popforecast==1) { free_ivector(popage,0,AGESUP); @@ -3187,10 +3220,10 @@ while((c=getc(ficpar))=='#' && c!= EOF){ covar=matrix(0,NCOVMAX,1,n); - cptcovn=0; + cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/ if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; - ncovmodel=2+cptcovn; + ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ /* Read guess parameters */ @@ -3228,7 +3261,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ fprintf(ficparo,"\n"); } - npar= (nlstate+ndeath-1)*nlstate*ncovmodel; + npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ p=param[1][1]; @@ -3387,12 +3420,12 @@ while((c=getc(ficpar))=='#' && c!= EOF){ Tvard=imatrix(1,15,1,2); Tage=ivector(1,15); - if (strlen(model) >1){ + if (strlen(model) >1){ /* If there is at least 1 covariate */ j=0, j1=0, k1=1, k2=1; - j=nbocc(model,'+'); - j1=nbocc(model,'*'); - cptcovn=j+1; - cptcovprod=j1; + j=nbocc(model,'+'); /* j=Number of '+' */ + j1=nbocc(model,'*'); /* j1=Number of '*' */ + cptcovn=j+1; + cptcovprod=j1; /*Number of products */ strcpy(modelsav,model); if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){ @@ -3401,6 +3434,8 @@ while((c=getc(ficpar))=='#' && c!= EOF){ goto end; } + /* This loop fill the array Tvar from the string 'model'.*/ + for(i=(j+1); i>=1;i--){ cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */ if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyze it */ @@ -3450,11 +3485,15 @@ while((c=getc(ficpar))=='#' && c!= EOF){ } /* end of loop + */ } /* end model */ + /*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products. + If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/ + /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]); printf("cptcovprod=%d ", cptcovprod); fprintf(ficlog,"cptcovprod=%d ", cptcovprod); - scanf("%d ",i);*/ - fclose(fic); + + scanf("%d ",i); + fclose(fic);*/ /* if(mle==1){*/ if (weightopt != 1) { /* Maximisation without weights*/ @@ -3525,8 +3564,8 @@ while((c=getc(ficpar))=='#' && c!= EOF){ } } -printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); - fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); + printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); + fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); free_vector(severity,1,maxwav); free_imatrix(outcome,1,maxwav+1,1,n); @@ -3545,13 +3584,15 @@ printf("Total number of individuals= %d, /* Concatenates waves */ concatwav(wav, dh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); + /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ Tcode=ivector(1,100); nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); ncodemax[1]=1; if (cptcovn > 0) tricode(Tvar,nbcode,imx); - codtab=imatrix(1,100,1,10); + codtab=imatrix(1,100,1,10); /* Cross tabulation to get the order of + the estimations*/ h=0; m=pow(2,cptcoveff); @@ -3704,8 +3745,8 @@ printf("Total number of individuals= %d, ungetc(c,ficpar); fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav); - fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,&mobilav); - fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,&mobilav); + fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); + fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); while((c=getc(ficpar))=='#' && c!= EOF){ ungetc(c,ficpar); @@ -3748,14 +3789,20 @@ while((c=getc(ficpar))=='#' && c!= EOF){ fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); - freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); + /*------------ gnuplot -------------*/ - strcpy(optionfilegnuplot,optionfilefiname); - strcat(optionfilegnuplot,".gp"); - if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { - printf("Problem with file %s",optionfilegnuplot); - } - fclose(ficgp); + strcpy(optionfilegnuplot,optionfilefiname); + strcat(optionfilegnuplot,".gp"); + if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { + printf("Problem with file %s",optionfilegnuplot); + } + else{ + fprintf(ficgp,"\n# %s\n", version); + fprintf(ficgp,"# %s\n", optionfilegnuplot); + fprintf(ficgp,"set missing 'NaNq'\n"); +} + fclose(ficgp); printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p); /*--------- index.htm --------*/ @@ -3792,17 +3839,17 @@ Interval (in months) between two waves: fclose(ficres); - /*--------------- Prevalence limit --------------*/ + /*--------------- Prevalence limit (stable prevalence) --------------*/ strcpy(filerespl,"pl"); strcat(filerespl,fileres); if((ficrespl=fopen(filerespl,"w"))==NULL) { - printf("Problem with Prev limit resultfile: %s\n", filerespl);goto end; - fprintf(ficlog,"Problem with Prev limit resultfile: %s\n", filerespl);goto end; + printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end; + fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end; } - printf("Computing prevalence limit: result on file '%s' \n", filerespl); - fprintf(ficlog,"Computing prevalence limit: result on file '%s' \n", filerespl); - fprintf(ficrespl,"#Prevalence limit\n"); + printf("Computing stable prevalence: result on file '%s' \n", filerespl); + fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl); + fprintf(ficrespl,"#Stable prevalence \n"); fprintf(ficrespl,"#Age "); for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); fprintf(ficrespl,"\n"); @@ -3948,11 +3995,17 @@ Interval (in months) between two waves: } printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); + calagedate=-1; + prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); - if (mobilav==1) { + + if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - movingaverage(probs, ageminpar, fage, mobaverage); + if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ + fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); + printf(" Error in movingaverage mobilav=%d\n",mobilav); + } } k=0; @@ -3994,10 +4047,10 @@ Interval (in months) between two waves: for(age=bage; age <=fage ;age++){ prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); if (popbased==1) { - if(mobilav !=1){ + if(mobilav ==0){ for(i=1; i<=nlstate;i++) prlim[i][i]=probs[(int)age][i][k]; - }else{ /* mobilav=1 */ + }else{ /* mobilav */ for(i=1; i<=nlstate;i++) prlim[i][i]=mobaverage[(int)age][i][k]; } @@ -4032,15 +4085,15 @@ free_matrix(mint,1,maxwav,1,n); fclose(ficpar); free_vector(epj,1,nlstate+1); - /*------- Variance limit prevalence------*/ + /*------- Variance of stable prevalence------*/ strcpy(fileresvpl,"vpl"); strcat(fileresvpl,fileres); if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { - printf("Problem with variance prev lim resultfile: %s\n", fileresvpl); + printf("Problem with variance of stable prevalence resultfile: %s\n", fileresvpl); exit(0); } - printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl); + printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl); k=0; for(cptcov=1;cptcov<=i1;cptcov++){ @@ -4075,7 +4128,7 @@ free_matrix(mint,1,maxwav,1,n); free_vector(delti,1,npar); free_matrix(agev,1,maxwav,1,imx); free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); - if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); fprintf(fichtm,"\n"); fclose(fichtm); @@ -4109,9 +4162,10 @@ free_matrix(mint,1,maxwav,1,n); strcpy(plotcmd,GNUPLOTPROGRAM); strcat(plotcmd," "); strcat(plotcmd,optionfilegnuplot); + printf("Starting: %s\n",plotcmd);fflush(stdout); system(plotcmd); -#ifdef windows + /*#ifdef windows*/ while (z[0] != 'q') { /* chdir(path); */ printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: "); @@ -4121,7 +4175,7 @@ free_matrix(mint,1,maxwav,1,n); else if (z[0] == 'g') system(plotcmd); else if (z[0] == 'q') exit(0); } -#endif + /*#endif */ }