--- imach/src/imach.c 2003/05/16 16:49:47 1.77 +++ imach/src/imach.c 2003/06/13 21:44:43 1.84 @@ -1,4 +1,21 @@ -/* $Id: imach.c,v 1.77 2003/05/16 16:49:47 brouard Exp $ +/* $Id: imach.c,v 1.84 2003/06/13 21:44:43 brouard Exp $ + $State: Exp $ + $Log: imach.c,v $ + Revision 1.84 2003/06/13 21:44:43 brouard + * imach.c (Repository): Replace "freqsummary" at a correct + place. It differs from routine "prevalence" which may be called + many times. Probs is memory consuming and must be used with + parcimony. + Version 0.95a2 (should output exactly the same maximization than 0.8a2) + + Revision 1.83 2003/06/10 13:39:11 lievre + *** empty log message *** + + Revision 1.82 2003/06/05 15:57:20 brouard + Add log in imach.c and fullversion number is now printed. + +*/ +/* Interpolated Markov Chain Short summary of the programme: @@ -58,6 +75,7 @@ read parameterfile read datafile concatwav + freqsummary if (mle >= 1) mlikeli print results files @@ -119,7 +137,11 @@ #define ODIRSEPARATOR '\\' #endif -char version[80]="Imach version 0.95a, May 2003, INED-EUROREVES "; +/* $Id: imach.c,v 1.84 2003/06/13 21:44:43 brouard Exp $ */ +/* $State: Exp $ */ + +char version[]="Imach version 0.95a2, June 2003, INED-EUROREVES "; +char fullversion[]="$Revision: 1.84 $ $Date: 2003/06/13 21:44:43 $"; int erreur; /* Error number */ int nvar; int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; @@ -1163,7 +1185,42 @@ double func( double *x) ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; } /* end of wave */ } /* end of individual */ - }else{ /* ml=4 no inter-extrapolation */ + }else if (mle==4){ /* ml=4 no inter-extrapolation */ + for (i=1,ipmx=0, sw=0.; i<=imx; i++){ + for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; + for(mi=1; mi<= wav[i]-1; mi++){ + for (ii=1;ii<=nlstate+ndeath;ii++) + for (j=1;j<=nlstate+ndeath;j++){ + oldm[ii][j]=(ii==j ? 1.0 : 0.0); + savm[ii][j]=(ii==j ? 1.0 : 0.0); + } + for(d=0; d nlstate){ + lli=log(out[s1][s2] - savm[s1][s2]); + }else{ + lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ + } + ipmx +=1; + sw += weight[i]; + ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; + /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/ + } /* end of wave */ + } /* end of individual */ + }else{ /* ml=5 no inter-extrapolation no jackson =0.8a */ for (i=1,ipmx=0, sw=0.; i<=imx; i++){ for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; for(mi=1; mi<= wav[i]-1; mi++){ @@ -1185,16 +1242,20 @@ double func( double *x) oldm=newm; } /* end mult */ + s1=s[mw[mi][i]][i]; + s2=s[mw[mi+1][i]][i]; lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ ipmx +=1; sw += weight[i]; ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; + /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/ } /* end of wave */ } /* end of individual */ } /* End of if */ for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; /* printf("l1=%f l2=%f ",ll[1],ll[2]); */ l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */ + /*exit(0); */ return -l; } @@ -1490,7 +1551,7 @@ void lubksb(double **a, int n, int *indx } /************ Frequencies ********************/ -void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2) +void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint) { /* Some frequencies */ int i, m, jk, k1,i1, j1, bool, z1,z2,j; @@ -1544,7 +1605,7 @@ void freqsummary(char fileres[], int ia if (bool==1){ for(m=firstpass; m<=lastpass; m++){ k2=anint[m][i]+(mint[m][i]/12.); - if ((k2>=dateprev1) && (k2<=dateprev2)) { + /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ if(agev[m][i]==0) agev[m][i]=iagemax+1; if(agev[m][i]==1) agev[m][i]=iagemax+2; if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i]; @@ -1557,12 +1618,12 @@ void freqsummary(char fileres[], int ia dateintsum=dateintsum+k2; k2cpt++; } - } + /*}*/ } } } - fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); + /* 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 "); @@ -1623,7 +1684,7 @@ void freqsummary(char fileres[], int ia if( i <= iagemax){ if(pos>=1.e-5){ fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); - probs[i][jk][j1]= pp[jk]/pos; + /*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 @@ -1656,7 +1717,7 @@ void freqsummary(char fileres[], int ia } /************ Prevalence ********************/ -void prevalence(double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) +void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) { /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people in each health status at the date of interview (if between dateprev1 and dateprev2). @@ -1799,7 +1860,7 @@ void concatwav(int wav[], int **dh, int sum=sum+j; /*if (j<0) printf("j=%d num=%d \n",j,i);*/ /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ - if(j<0)printf("Warning! Negative delay (%d) between waves %d and %d of individual at line %d who is aged %.1f with statuses %d %d\n ",j,mw[mi][i],mw[mi+1][i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + if(j<0)printf("Error! Negative delay (%d to death) between waves %d and %d of individual %d at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); } } else{ @@ -1810,7 +1871,7 @@ void concatwav(int wav[], int **dh, int 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]);*/ - if(j<0)printf("Warning! Negative delay (%d to death) between waves %d and %d of individual at line %d who is aged %.1f with statuses %d %d\n ",j,mw[mi][i],mw[mi+1][i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + if(j<0)printf("Error! Negative delay (%d) between waves %d and %d of individual %d at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); sum=sum+j; } jk= j/stepm; @@ -2390,7 +2451,7 @@ void varevsij(char optionfilefiname[], d fclose(ficresprobmorprev); fclose(ficgp); fclose(fichtm); -} +} /* end varevsij */ /************ Variance of prevlim ******************/ void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij) @@ -2987,7 +3048,7 @@ m=pow(2,cptcoveff); fprintf(ficgp,"\nset out \"p%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1); - for (i=1; i<= nlstate ; i ++) + for (i=1; i< nlstate ; i ++) fprintf(ficgp,"+$%d",k+i+1); fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1); @@ -3126,7 +3187,7 @@ prevforecast(char fileres[], double anpr char fileresf[FILENAMELENGTH]; agelim=AGESUP; - prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); + prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); strcpy(fileresf,"f"); strcat(fileresf,fileres); @@ -3249,7 +3310,7 @@ populforecast(char fileres[], double anp agelim=AGESUP; calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM; - prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); + prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); strcpy(filerespop,"pop"); @@ -3391,7 +3452,7 @@ populforecast(char fileres[], double anp free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); fclose(ficrespop); -} +} /* End of popforecast */ /***********************************************/ /**************** Main Program *****************/ @@ -3454,7 +3515,7 @@ int main(int argc, char *argv[]) gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ getcwd(pathcd, size); - printf("\n%s",version); + printf("\n%s\n%s",version,fullversion); if(argc <=1){ printf("\nEnter the parameter file name: "); scanf("%s",pathtot); @@ -3832,7 +3893,7 @@ int main(int argc, char *argv[]) if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){ printf("Error! Month of death of individual %d on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); fprintf(ficlog,"Error! Month of death of individual %d on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); - s[m][i]=-1; + s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */ } } } @@ -3954,6 +4015,7 @@ int main(int argc, char *argv[]) /* Calculates basic frequencies. Computes observed prevalence at single age and prints on file fileres'p'. */ + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ @@ -3997,7 +4059,7 @@ int main(int argc, char *argv[]) } } } - if(mle==1){ + if(mle!=0){ /* Computing hessian and covariance matrix */ ftolhess=ftol; /* Usually correct */ hesscov(matcov, p, npar, delti, ftolhess, func); @@ -4128,8 +4190,8 @@ int main(int argc, char *argv[]) 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); - probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); - freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); + /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ /*------------ gnuplot -------------*/ strcpy(optionfilegnuplot,optionfilefiname); @@ -4293,6 +4355,7 @@ Interval (in months) between two waves: fclose(ficrespij); + probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); /*---------- Forecasting ------------------*/ /*if((stepm == 1) && (strcmp(model,".")==0)){*/ @@ -4340,7 +4403,7 @@ Interval (in months) between two waves: fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */ - prevalence(agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); + prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); /* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\ ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass); */