--- imach/src/imach.c 2002/03/13 17:19:16 1.34 +++ imach/src/imach.c 2002/03/26 17:08:39 1.35 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.34 2002/03/13 17:19:16 brouard Exp $ +/* $Id: imach.c,v 1.35 2002/03/26 17:08:39 lievre Exp $ Interpolated Markov Chain Short summary of the programme: @@ -56,7 +56,8 @@ #include #define MAXLINE 256 -#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot" +#define GNUPLOTPROGRAM "wgnuplot" +/*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ #define FILENAMELENGTH 80 /*#define DEBUG*/ #define windows @@ -686,16 +687,15 @@ double **prevalim(double **prlim, int nl for (k=1; k<=cptcovn;k++) { cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; - /*printf("ij=%d Tvar[k]=%d nbcode=%d cov=%lf\n",ij, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k]);*/ + /* printf("ij=%d k=%d Tvar[k]=%d nbcode=%d cov=%lf codtab[ij][Tvar[k]]=%d \n",ij,k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], codtab[ij][Tvar[k]]);*/ } - for (k=1; k<=cptcovage;k++) - cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; + for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; for (k=1; k<=cptcovprod;k++) cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/ - + /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/ out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); savm=oldm; @@ -1178,14 +1178,14 @@ 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,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; double ***freq; /* Frequencies */ double *pp; double pos, k2, dateintsum=0,k2cpt=0; FILE *ficresp; char fileresp[FILENAMELENGTH]; - + pp=vector(1,nlstate); probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); strcpy(fileresp,"p"); @@ -1196,113 +1196,116 @@ void freqsummary(char fileres[], int ag } freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3); j1=0; - + j=cptcoveff; if (cptcovn<1) {j=1;ncodemax[1]=1;} - + for(k1=1; k1<=j;k1++){ - for(i1=1; i1<=ncodemax[k1];i1++){ - j1++; - /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); - scanf("%d", i);*/ - for (i=-1; i<=nlstate+ndeath; i++) - for (jk=-1; jk<=nlstate+ndeath; jk++) - for(m=agemin; m <= agemax+3; m++) - freq[i][jk][m]=0; - - dateintsum=0; - k2cpt=0; - for (i=1; i<=imx; i++) { - bool=1; - if (cptcovn>0) { - for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) - bool=0; - } - if (bool==1) { - for(m=firstpass; m<=lastpass; m++){ - k2=anint[m][i]+(mint[m][i]/12.); - if ((k2>=dateprev1) && (k2<=dateprev2)) { - 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]; - freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i]; - if ((agev[m][i]>1) && (agev[m][i]< (agemax+3))) { - dateintsum=dateintsum+k2; - k2cpt++; - } - - } - } - } - } + for(i1=1; i1<=ncodemax[k1];i1++){ + j1++; + /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); + scanf("%d", i);*/ + for (i=-1; i<=nlstate+ndeath; i++) + for (jk=-1; jk<=nlstate+ndeath; jk++) + for(m=agemin; m <= agemax+3; m++) + freq[i][jk][m]=0; + + dateintsum=0; + k2cpt=0; + for (i=1; i<=imx; i++) { + bool=1; + if (cptcovn>0) { + for (z1=1; z1<=cptcoveff; z1++) + if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) + bool=0; + } + if (bool==1) { + for(m=firstpass; m<=lastpass; m++){ + k2=anint[m][i]+(mint[m][i]/12.); + if ((k2>=dateprev1) && (k2<=dateprev2)) { + if(agev[m][i]==0) agev[m][i]=agemax+1; + if(agev[m][i]==1) agev[m][i]=agemax+2; + if (m1) && (agev[m][i]< (agemax+3))) { + 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 "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - fprintf(ficresp, "**********\n#"); - } - for(i=1; i<=nlstate;i++) - fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); - fprintf(ficresp, "\n"); - - for(i=(int)agemin; i <= (int)agemax+3; i++){ - if(i==(int)agemax+3) - printf("Total"); - else - 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]; - } - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pos=0; m <=0 ; m++) - pos += freq[jk][m][i]; - if(pp[jk]>=1.e-10) - printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); - else - printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); - } + 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]]); + fprintf(ficresp, "**********\n#"); + } + for(i=1; i<=nlstate;i++) + fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); + fprintf(ficresp, "\n"); + + for(i=(int)agemin; i <= (int)agemax+3; i++){ + if(i==(int)agemax+3) + printf("Total"); + else + 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]; + } + for(jk=1; jk <=nlstate ; jk++){ + for(m=-1, pos=0; m <=0 ; m++) + pos += freq[jk][m][i]; + if(pp[jk]>=1.e-10) + printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + else + printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); + } - 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; 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++){ - if(pos>=1.e-5) - printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); - else - printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); - if( i <= (int) agemax){ - if(pos>=1.e-5){ - fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, 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]);*/ + for(jk=1,pos=0; jk <=nlstate ; jk++) + pos += pp[jk]; + for(jk=1; jk <=nlstate ; jk++){ + if(pos>=1.e-5) + printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); + else + printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); + if( i <= (int) agemax){ + if(pos>=1.e-5){ + fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, 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 + fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos); + } } - else - fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos); + + for(jk=-1; jk <=nlstate+ndeath; jk++) + for(m=-1; m <=nlstate+ndeath; m++) + if(freq[jk][m][i] !=0 ) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); + if(i <= (int) agemax) + fprintf(ficresp,"\n"); + printf("\n"); } } - for(jk=-1; jk <=nlstate+ndeath; jk++) - for(m=-1; m <=nlstate+ndeath; m++) - if(freq[jk][m][i] !=0 ) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); - if(i <= (int) agemax) - fprintf(ficresp,"\n"); - printf("\n"); - } - } - } + } dateintmean=dateintsum/k2cpt; fclose(ficresp); free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); free_vector(pp,1,nlstate); - + /* End of Freq */ } @@ -1346,7 +1349,7 @@ void prevalence(int agemin, float agemax if ((k2>=dateprev1) && (k2<=dateprev2)) { 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]+1-((int)calagedate %12)/12.)] += weight[i]; + if (m ncodemax[j]) break; @@ -1528,7 +1532,7 @@ void evsij(char fileres[], double ***eij { /* Health expectancies */ int i, j, nhstepm, hstepm, h, nstepm, k; - double age, agelim,hf; + double age, agelim, hf; double ***p3mat; fprintf(ficreseij,"# Health expectancies\n"); @@ -1589,13 +1593,13 @@ void varevsij(char fileres[], double *** /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ double **newm; double **dnewm,**doldm; - int i, j, nhstepm, hstepm, h; + int i, j, nhstepm, hstepm, h, nstepm, kk; int k, cptcode; double *xp; double **gp, **gm; double ***gradg, ***trgradg; double ***p3mat; - double age,agelim; + double age,agelim, hf; int theta; fprintf(ficresvij,"# Covariances of life expectancies\n"); @@ -1609,13 +1613,25 @@ void varevsij(char fileres[], double *** dnewm=matrix(1,nlstate,1,npar); doldm=matrix(1,nlstate,1,nlstate); - hstepm=1*YEARM; /* Every year of age */ - hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */ + kk=1; /* For example stepm=6 months */ + hstepm=kk*YEARM; /* (a) Every k years of age (in months), for example every k=2 years 24 m */ + hstepm=stepm; /* or (b) We decided to compute the life expectancy with the smallest unit */ + /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. + nhstepm is the number of hstepm from age to agelim + nstepm is the number of stepm from age to agelin. + Look at hpijx to understand the reason of that which relies in memory size + and note for a fixed period like k years */ + /* We decided (b) to get a life expectancy respecting the most precise curvature of the + survival function given by stepm (the optimization length). Unfortunately it + means that if the survival funtion is printed only each two years of age and if + you sum them up and add 1 year (area under the trapezoids) you won't get the same + results. So we changed our mind and took the option of the best precision. + */ + hstepm=hstepm/stepm; /* Typically in stepm units, if k= 2 years, = 2/6 months = 4 */ agelim = AGESUP; for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ - nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ - if (stepm >= YEARM) hstepm=1; - nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ + nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ + nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); gradg=ma3x(0,nhstepm,1,npar,1,nlstate); gp=matrix(0,nhstepm,1,nlstate); @@ -1670,24 +1686,25 @@ void varevsij(char fileres[], double *** for(theta=1; theta <=npar; theta++) trgradg[h][j][theta]=gradg[h][theta][j]; + hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ for(i=1;i<=nlstate;i++) for(j=1;j<=nlstate;j++) vareij[i][j][(int)age] =0.; + for(h=0;h<=nhstepm;h++){ for(k=0;k<=nhstepm;k++){ matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov); matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]); for(i=1;i<=nlstate;i++) for(j=1;j<=nlstate;j++) - vareij[i][j][(int)age] += doldm[i][j]; + vareij[i][j][(int)age] += doldm[i][j]*hf*hf; } } - h=1; - if (stepm >= YEARM) h=stepm/YEARM; + fprintf(ficresvij,"%.0f ",age ); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++){ - fprintf(ficresvij," %.4f", h*vareij[i][j][(int)age]); + fprintf(ficresvij," %.4f", vareij[i][j][(int)age]); } fprintf(ficresvij,"\n"); free_matrix(gp,0,nhstepm,1,nlstate); @@ -1892,7 +1909,11 @@ fclose(ficresprob); } /******************* Printing html file ***********/ -void printinghtml(char fileres[], char title[], char datafile[], int firstpass, int lastpass, int stepm, int weightopt, char model[],int imx,int jmin, int jmax, double jmeanint,char optionfile[],char optionfilehtm[],char rfileres[] ){ +void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \ + int lastpass, int stepm, int weightopt, char model[],\ + int imx,int jmin, int jmax, double jmeanint,char optionfile[], \ + char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\ + char version[], int popforecast ){ int jj1, k1, i1, cpt; FILE *fichtm; /*char optionfilehtm[FILENAMELENGTH];*/ @@ -1903,26 +1924,32 @@ void printinghtml(char fileres[], char t printf("Problem with %s \n",optionfilehtm), exit(0); } - fprintf(fichtm,"
    Imach, Version 0.8
    -Title=%s
    Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
    - -Total number of observations=%d
    -Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
    + fprintf(fichtm," Imach, Version %s
    \n +Title=%s
    Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
    \n +\n +Total number of observations=%d
    \n +Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
    \n
    -
  • Outputs files

    \n - - Observed prevalence in each state: p%s
    \n -- Estimated parameters and the covariance matrix: %s
    - - Stationary prevalence in each state: pl%s
    - - Transition probabilities: pij%s
    - - Copy of the parameter file: o%s
    - - 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
    - - Prevalences forecasting: f%s
    - - Population forecasting (if popforecast=1): pop%s
    -
    ",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres); - +
    • Outputs files
      \n + - Copy of the parameter file: o%s
      \n + - Gnuplot file name: %s
      \n + - Observed prevalence in each state: p%s
      \n + - Stationary prevalence in each state: pl%s
      \n + - Transition probabilities: pij%s
      \n + - Life expectancies by age and initial health status: e%s
      \n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres); + + fprintf(fichtm,"\n + - Parameter file with estimated parameters and the covariance matrix: %s
      \n + - Variances of life expectancies by age and initial health status: v%s
      \n + - Health expectancies with their variances: t%s
      \n + - Standard deviation of stationary prevalences: vpl%s
      \n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres); + + if(popforecast==1) fprintf(fichtm,"\n + - Prevalences forecasting: f%s
      \n + - Population forecasting (if popforecast=1): pop%s
      \n +
      ",fileres,fileres,fileres,fileres); + else + fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)

      • \n",popforecast, stepm, model); fprintf(fichtm,"
    • Graphs

      • "); m=cptcoveff; @@ -1963,12 +1990,12 @@ fclose(fichtm); } /******************* Gnuplot file **************/ -void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double agemin, double agemaxpar, double fage , char pathc[], double p[]){ +void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ int m,cpt,k1,i,k,j,jk,k2,k3,ij,l; strcpy(optionfilegnuplot,optionfilefiname); - strcat(optionfilegnuplot,".plt"); + strcat(optionfilegnuplot,".gp.txt"); if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { printf("Problem with file %s",optionfilegnuplot); } @@ -1983,10 +2010,10 @@ m=pow(2,cptcoveff); for (k1=1; k1<= m ; k1 ++) { #ifdef windows - fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",agemin,fage,fileres,k1-1,k1-1); + fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1); #endif #ifdef unix -fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",agemin,fage,fileres); +fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres); #endif for (i=1; i<= nlstate ; i ++) { @@ -2013,7 +2040,7 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\" /*2 eme*/ for (k1=1; k1<= m ; k1 ++) { - fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",agemin,fage); + fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",ageminpar,fage); for (i=1; i<= nlstate+1 ; i ++) { k=2*i; @@ -2046,7 +2073,7 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\" for (k1=1; k1<= m ; k1 ++) { for (cpt=1; cpt<= nlstate ; cpt ++) { k=2+nlstate*(cpt-1); - fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",agemin,fage,fileres,k1-1,k1-1,k,cpt); + fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt); for (i=1; i< nlstate ; i ++) { fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+i,cpt,i+1); } @@ -2058,7 +2085,7 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\" for (k1=1; k1<= m ; k1 ++) { for (cpt=1; cpt=(agemin-((int)calagedate %12)/12.); agedeb--){ + for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){ nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); nhstepm = nhstepm/hstepm; @@ -2231,7 +2258,7 @@ calagedate=(anproj1+mproj1/12.+jproj1/36 for (h=0; h<=nhstepm; h++){ if (h==(int) (calagedate+YEARM*cpt)) { - fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm); + fprintf(ficresf,"\n %.f %.f ",anproj1+cpt,agedeb+h*hstepm/YEARM*stepm); } for(j=1; j<=nlstate+ndeath;j++) { kk1=0.;kk2=0; @@ -2260,7 +2287,7 @@ calagedate=(anproj1+mproj1/12.+jproj1/36 fclose(ficresf); } /************** Forecasting ******************/ -populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double agemin, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ +populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h; int *popage; @@ -2274,7 +2301,7 @@ populforecast(char fileres[], double anp agelim=AGESUP; calagedate=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM; - prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); + prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); strcpy(filerespop,"pop"); @@ -2288,7 +2315,7 @@ populforecast(char fileres[], double anp if (mobilav==1) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - movingaverage(agedeb, fage, agemin, mobaverage); + movingaverage(agedeb, fage, ageminpar, mobaverage); } stepsize=(int) (stepm+YEARM-1)/YEARM; @@ -2329,7 +2356,7 @@ populforecast(char fileres[], double anp for (cpt=0; cpt<=0;cpt++) { fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt); - for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ + for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){ nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); nhstepm = nhstepm/hstepm; @@ -2375,7 +2402,7 @@ populforecast(char fileres[], double anp for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) { fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt); - for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ + for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){ nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); nhstepm = nhstepm/hstepm; @@ -2421,7 +2448,7 @@ int main(int argc, char *argv[]) int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod; double agedeb, agefin,hf; - double agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; + double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; double fret; double **xi,tmp,delta; @@ -2466,7 +2493,7 @@ int main(int argc, char *argv[]) double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2; - char version[80]="Imach version 0.8, March 2002, INED-EUROREVES "; + char version[80]="Imach version 0.8a, March 2002, INED-EUROREVES "; char *alph[]={"a","a","b","c","d","e"}, str[4]; @@ -2479,7 +2506,7 @@ int main(int argc, char *argv[]) struct timeval start_time, end_time; gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ - + getcwd(pathcd, size); printf("\n%s",version); if(argc <=1){ @@ -2705,11 +2732,13 @@ while((c=getc(ficpar))=='#' && c!= EOF){ 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++) - if (covar[1][i]==0) 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++){ + if (s[4][i]==9) s[4][i]=-1; + 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); Tprod=ivector(1,15); @@ -2724,7 +2753,6 @@ while((c=getc(ficpar))=='#' && c!= EOF){ cptcovn=j+1; cptcovprod=j1; - strcpy(modelsav,model); if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){ printf("Error. Non available option model=%s ",model); @@ -2780,7 +2808,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ } } - /*printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]); + /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]); printf("cptcovprod=%d ", cptcovprod); scanf("%d ",i);*/ fclose(fic); @@ -2792,22 +2820,26 @@ while((c=getc(ficpar))=='#' && c!= EOF){ /*-calculation of age at interview from date of interview and age at death -*/ agev=matrix(1,maxwav,1,imx); - for (i=1; i<=imx; i++) - for(m=2; (m<= maxwav); m++) + for (i=1; i<=imx; i++) { + for(m=2; (m<= maxwav); m++) { if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){ anint[m][i]=9999; s[m][i]=-1; } - + if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1; + } + } + 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] == nlstate+1) { + if (s[m][i] >= nlstate+1) { if(agedc[i]>0) if(moisdc[i]!=99 && andc[i]!=9999) - agev[m][i]=agedc[i]; - else { + agev[m][i]=agedc[i]; + /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ + else { if (andc[i]!=9999){ printf("Warning negative age at death: %d line:%d\n",num[i],i); agev[m][i]=-1; @@ -2882,20 +2914,21 @@ printf("Total number of individuals= %d, for(j=1; j <= ncodemax[k]; j++){ for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){ h++; - if (h>m) h=1;codtab[h][k]=j; + if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j; + /* printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/ } } } } - - - /*for(i=1; i <=m ;i++){ - for(k=1; k <=cptcovn; k++){ - printf("i=%d k=%d %d %d",i,k,codtab[i][k], cptcoveff); - } - printf("\n"); - } - scanf("%d",i);*/ + /* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]); + codtab[1][2]=1;codtab[2][2]=2; */ + /* for(i=1; i <=m ;i++){ + for(k=1; k <=cptcovn; k++){ + printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); + } + printf("\n"); + } + scanf("%d",i);*/ /* Calculates basic frequencies. Computes observed prevalence at single age and prints on file fileres'p'. */ @@ -2989,16 +3022,16 @@ printf("Total number of individuals= %d, } ungetc(c,ficpar); - fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemaxpar, &bage, &fage); + fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&ageminpar,&agemaxpar, &bage, &fage); if (fage <= 2) { - bage = agemin; + bage = ageminpar; fage = agemaxpar; } fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); - fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemaxpar,bage,fage); - fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemaxpar,bage,fage); + fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",ageminpar,agemaxpar,bage,fage); + fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",ageminpar,agemaxpar,bage,fage); while((c=getc(ficpar))=='#' && c!= EOF){ ungetc(c,ficpar); @@ -3056,7 +3089,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ 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,agemaxpar,fage, pathc,p); + printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, ageminpar,agemaxpar,fage, pathc,p); /*------------ free_vector -------------*/ chdir(path); @@ -3072,7 +3105,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ /*--------- index.htm --------*/ - printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres); + printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast); /*--------------- Prevalence limit --------------*/ @@ -3095,7 +3128,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */ k=0; - agebase=agemin; + agebase=ageminpar; agelim=agemaxpar; ftolpl=1.e-10; i1=cptcoveff; @@ -3222,12 +3255,12 @@ while((c=getc(ficpar))=='#' && c!= EOF){ fprintf(ficreseij,"\n#****** "); for(j=1;j<=cptcoveff;j++) - fprintf(ficreseij,"V%d=%d ",j,nbcode[j][codtab[k][j]]); + fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); fprintf(ficreseij,"******\n"); fprintf(ficresvij,"\n#****** "); for(j=1;j<=cptcoveff;j++) - fprintf(ficresvij,"V%d=%d ",j,nbcode[j][codtab[k][j]]); + fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); fprintf(ficresvij,"******\n"); eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); @@ -3345,14 +3378,12 @@ while((c=getc(ficpar))=='#' && c!= EOF){ #ifdef windows while (z[0] != 'q') { - chdir(path); - printf("\nType e to edit output files, c to start again, and q for exiting: "); + /* chdir(path); */ + printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: "); scanf("%s",z); if (z[0] == 'c') system("./imach"); - else if (z[0] == 'e') { - chdir(path); - system(optionfilehtm); - } + else if (z[0] == 'e') system(optionfilehtm); + else if (z[0] == 'g') system(plotcmd); else if (z[0] == 'q') exit(0); } #endif