--- imach/src/imach.c 2003/01/28 17:23:35 1.66 +++ imach/src/imach.c 2003/02/04 12:40:59 1.68 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.66 2003/01/28 17:23:35 brouard Exp $ +/* $Id: imach.c,v 1.68 2003/02/04 12:40:59 lievre Exp $ Interpolated Markov Chain Short summary of the programme: @@ -919,7 +919,7 @@ double func( double *x) double sw; /* Sum of weights */ double lli; /* Individual log likelihood */ int s1, s2; - double bbh; + double bbh, survp; long ipmx; /*extern weight */ /* We are differentiating ll according to initial status */ @@ -972,6 +972,14 @@ double func( double *x) * is higher than the multiple of stepm and negative otherwise. */ /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ + /* if s2=-2 lli=out[1][1]+out[1][2];*/ + if (s2==-2) { + for (j=1,survp=0. ; j<=nlstate; j++) + survp += out[s1][j]; + lli= survp; + } + + else lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ /*if(lli ==000.0)*/ @@ -1416,7 +1424,7 @@ void freqsummary(char fileres[], int ag fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); exit(0); } - freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3); + freq= ma3x(-2,nlstate+ndeath,-2,nlstate+ndeath,agemin,agemax+3); j1=0; j=cptcoveff; @@ -1531,8 +1539,8 @@ void freqsummary(char fileres[], int ag } } - for(jk=-1; jk <=nlstate+ndeath; jk++) - for(m=-1; m <=nlstate+ndeath; m++) + for(jk=-2; jk <=nlstate+ndeath; jk++) + for(m=-2; m <=nlstate+ndeath; m++) if(freq[jk][m][i] !=0 ) { if(first==1) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); @@ -1549,7 +1557,7 @@ void freqsummary(char fileres[], int ag dateintmean=dateintsum/k2cpt; fclose(ficresp); - free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); + free_ma3x(freq,-2,nlstate+ndeath,-2,nlstate+ndeath,(int) agemin,(int) agemax+3); free_vector(pp,1,nlstate); /* End of Freq */ @@ -1566,7 +1574,7 @@ void prevalence(int agemin, float agemax pp=vector(1,nlstate); - freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3); + freq=ma3x(-2,nlstate+ndeath,-2,nlstate+ndeath,agemin,agemax+3); j1=0; j=cptcoveff; @@ -1576,8 +1584,8 @@ 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 (i=-2; i<=nlstate+ndeath; i++) + for (jk=-2; jk<=nlstate+ndeath; jk++) for(m=agemin; m <= agemax+3; m++) freq[i][jk][m]=0; @@ -1634,7 +1642,7 @@ void prevalence(int agemin, float agemax } /* end k1 */ - free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); + free_ma3x(freq,-2,nlstate+ndeath,-2,nlstate+ndeath,(int) agemin,(int) agemax+3); free_vector(pp,1,nlstate); } /* End of Freq */ @@ -1664,7 +1672,7 @@ void concatwav(int wav[], int **dh, int mi=0; m=firstpass; while(s[m][i] <= nlstate){ - if(s[m][i]>=1) + if(s[m][i]>=1 || s[m][i]==-2) mw[++mi][i]=m; if(m >=lastpass) break; @@ -1704,10 +1712,12 @@ void concatwav(int wav[], int **dh, int if (j <= jmin) jmin=j; 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);*/ } } else{ j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); + /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ k=k+1; if (j >= jmax) jmax=j; else if (j <= jmin)jmin=j; @@ -1742,7 +1752,6 @@ void concatwav(int wav[], int **dh, int bh[mi][i]=ju; /* At least one step */ printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i); } - if(i==298 || i==287 || i==763 ||i==1061)printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d",bh[mi][i],ju,jl,dh[mi][i],jk,stepm); } } /* end if mle */ } /* end wave */ @@ -2058,8 +2067,8 @@ void varevsij(char optionfilefiname[], d exit(0); } else{ - fprintf(fichtm,"\n

  • Computing probabilities of dying as a weighted average (i.e global mortality independent of initial healh state)

    • \n"); - fprintf(fichtm,"\n
    %s (à revoir)
    \n",digitp); + fprintf(fichtm,"\n

  • Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)

    • \n"); + fprintf(fichtm,"\n
    %s
    \n",digitp); } varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); @@ -2136,7 +2145,7 @@ void varevsij(char optionfilefiname[], d as a weighted average of prlim. */ for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){ - for(i=1; i<= nlstate; i++) + for(i=1,gpp[j]=0.; i<= nlstate; i++) gpp[j] += prlim[i][i]*p3mat[i][j][1]; } /* end probability of death */ @@ -2167,8 +2176,8 @@ void varevsij(char optionfilefiname[], d as a weighted average of prlim. */ for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ - for(i=1; i<= nlstate; i++) - gmp[j] += prlim[i][i]*p3mat[i][j][1]; + for(i=1,gmp[j]=0.; i<= nlstate; i++) + gmp[j] += prlim[i][i]*p3mat[i][j][1]; } /* end probability of death */ @@ -2176,6 +2185,7 @@ void varevsij(char optionfilefiname[], d for(h=0; h<=nhstepm; h++){ gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; } + for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */ gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta]; } @@ -2190,8 +2200,9 @@ void varevsij(char optionfilefiname[], d trgradg[h][j][theta]=gradg[h][theta][j]; for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */ - for(theta=1; theta <=npar; theta++) + for(theta=1; theta <=npar; theta++) { trgradgp[j][theta]=gradgp[theta][j]; + } hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ for(i=1;i<=nlstate;i++) @@ -2211,9 +2222,12 @@ void varevsij(char optionfilefiname[], d /* pptj */ matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov); matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp); - for(j=nlstate+1;j<=nlstate+ndeath;j++) - for(i=nlstate+1;i<=nlstate+ndeath;i++) + + for(j=nlstate+1;j<=nlstate+ndeath;j++) + for(i=nlstate+1;i<=nlstate+ndeath;i++){ varppt[j][i]=doldmp[j][i]; + } + /* end ppptj */ /* x centered again */ hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); @@ -2233,8 +2247,8 @@ void varevsij(char optionfilefiname[], d computed over hstepm (estepm) matrices product = hstepm*stepm months) as a weighted average of prlim. */ - for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ - for(i=1; i<= nlstate; i++) + for(j=nlstate+1;j<=nlstate+ndeath;j++){ + for(i=1,gmp[j]=0.;i<= nlstate; i++) gmp[j] += prlim[i][i]*p3mat[i][j][1]; } /* end probability of death */ @@ -2267,9 +2281,12 @@ void varevsij(char optionfilefiname[], d fprintf(ficgp,"\nset noparametric;set nolabel; set ter png small;set size 0.65, 0.65"); /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */ fprintf(ficgp,"\n set log y; set nolog x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); - fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); - fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); - fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); +/* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ +/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ +/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ + fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l 1 ",fileresprobmorprev); + fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95\%% interval\" w l 2 ",fileresprobmorprev); + fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l 2 ",fileresprobmorprev); fprintf(fichtm,"\n
    File (multiple files are possible if covariates are present): %s\n",fileresprobmorprev,fileresprobmorprev); fprintf(fichtm,"\n
    Probability is computed over estepm=%d months.

    \n", estepm,digitp,digit); /* fprintf(fichtm,"\n
    Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year

    \n", stepm,YEARM,digitp,digit); @@ -3701,11 +3718,10 @@ int main(int argc, char *argv[]) for (i=1; i<=imx; i++) { agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); for(m=1; (m<= maxwav); m++){ - if(s[m][i] >0){ + if(s[m][i] >0 || s[m][i]==-2){ if (s[m][i] >= nlstate+1) { if(agedc[i]>0) - if(moisdc[i]!=99 && andc[i]!=9999) - agev[m][i]=agedc[i]; + if(moisdc[i]!=99 && andc[i]!=9999) agev[m][i]=agedc[i]; /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ else { if (andc[i]!=9999){ @@ -3767,7 +3783,8 @@ int main(int argc, char *argv[]) dh=imatrix(1,lastpass-firstpass+1,1,imx); bh=imatrix(1,lastpass-firstpass+1,1,imx); mw=imatrix(1,lastpass-firstpass+1,1,imx); - + + /* Concatenates waves */ concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); @@ -3976,6 +3993,7 @@ 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); + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); /*------------ gnuplot -------------*/ @@ -4138,7 +4156,7 @@ Interval (in months) between two waves: /*---------- Forecasting ------------------*/ - if((stepm == 1) && (strcmp(model,".")==0)){ + if((stepm == 1) && (strcmp(model,".")==0)){ prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1); if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1); } @@ -4146,7 +4164,7 @@ Interval (in months) between two waves: erreur=108; printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); - } + } /*---------- Health expectancies and variances ------------*/ @@ -4170,6 +4188,7 @@ Interval (in months) between two waves: printf("Computing Health Expectancies: result on file '%s' \n", filerese); fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese); + strcpy(fileresv,"v"); strcat(fileresv,fileres); if((ficresvij=fopen(fileresv,"w"))==NULL) {