--- imach096d/src/imach.c 2001/05/09 14:25:42 1.10 +++ imach096d/src/imach.c 2002/02/21 18:42:24 1.21 @@ -8,7 +8,7 @@ Health expectancies are computed from the transistions observed between waves and are computed for each degree of severity of disability (number of life states). More degrees you consider, more time is necessary to - reach the Maximum Likelihood of the parameters involved in the model. + reach the Maximum Likelihood of the parameters involved in the model. The simplest model is the multinomial logistic model where pij is the probabibility to be observed in state j at the second wave conditional to be observed in state i at the first wave. Therefore the model is: @@ -22,7 +22,7 @@ delay between waves is not identical for each individual, or if some individual missed an interview, the information is not rounded or lost, but taken into account using an interpolation or extrapolation. - hPijx is the probability to be + hPijx is the probability to be observed in state i at age x+h conditional to the observed state i at age x. The delay 'h' can be split into an exact number (nh*stepm) of unobserved intermediate states. This elementary transition (by month or @@ -67,12 +67,14 @@ #define AGEBASE 40 +int erreur; /* Error number */ int nvar; int cptcovn, cptcovage=0, cptcoveff=0,cptcov; int npar=NPARMAX; int nlstate=2; /* Number of live states */ int ndeath=1; /* Number of dead states */ int ncovmodel, ncov; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */ +int popbased=0; int *wav; /* Number of waves for this individuual 0 is possible */ int maxwav; /* Maxim number of waves */ @@ -83,8 +85,8 @@ 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; -FILE *ficgp, *fichtm; +FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf; +FILE *ficgp, *fichtm,*ficresprob,*ficpop; FILE *ficreseij; char filerese[FILENAMELENGTH]; FILE *ficresvij; @@ -127,7 +129,8 @@ int stepm; int m,nb; int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage; double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint; -double **pmmij; +double **pmmij, ***probs, ***mobaverage; +double dateintmean=0; double *weight; int **s; /* Status */ @@ -693,7 +696,7 @@ double **prevalim(double **prlim, int nl } } -/*************** transition probabilities **********/ +/*************** transition probabilities ***************/ double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate ) { @@ -716,9 +719,11 @@ double **pmij(double **ps, double *cov, s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc]; /*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/ } - ps[i][j]=s2; + ps[i][j]=(s2); } } + /*ps[3][2]=1;*/ + for(i=1; i<= nlstate; i++){ s1=0; for(j=1; ji) { printf(".%d%d",i,j);fflush(stdout); hess[i][j]=hessij(p,delti,i,j); - hess[j][i]=hess[i][j]; + hess[j][i]=hess[i][j]; + /*printf(" %lf ",hess[i][j]);*/ } } } @@ -1032,7 +1039,7 @@ double hessii( double x[], double delta, } } delti[theta]=delts; - return res; + return res; } @@ -1145,18 +1152,18 @@ 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) +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) { /* Some frequencies */ - int i, m, jk, k1, i1, j1, bool, z1,z2,j; + int i, m, jk, k1,i1, j1, bool, z1,z2,j; double ***freq; /* Frequencies */ double *pp; - double pos; + 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"); strcat(fileresp,fileres); if((ficresp=fopen(fileresp,"w"))==NULL) { @@ -1178,7 +1185,9 @@ void freqsummary(char fileres[], int ag 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) { @@ -1186,12 +1195,20 @@ void freqsummary(char fileres[], int ag if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) bool=0; } - if (bool==1) { - for(m=firstpass; m<=lastpass-1; m++){ - 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 (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++; + } + + } } } } @@ -1211,7 +1228,7 @@ void freqsummary(char fileres[], int ag 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]; + pp[jk] += freq[jk][m][i]; } for(jk=1; jk <=nlstate ; jk++){ for(m=-1, pos=0; m <=0 ; m++) @@ -1221,10 +1238,12 @@ void freqsummary(char fileres[], int ag else printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); } - for(jk=1; jk <=nlstate ; jk++){ - for(m=1, pp[jk]=0; m <=nlstate+ndeath; m++) + + 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++){ @@ -1233,8 +1252,11 @@ void freqsummary(char fileres[], int ag else printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); if( i <= (int) agemax){ - if(pos>=1.e-5) + 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); } @@ -1248,11 +1270,95 @@ 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_vector(pp,1,nlstate); + /* End of Freq */ +} + +/************ Prevalence ********************/ +void prevalence(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 calagedate) +{ /* Some frequencies */ + + int i, m, jk, k1, i1, j1, bool, z1,z2,j; + double ***freq; /* Frequencies */ + double *pp; + double pos, k2; + + pp=vector(1,nlstate); + probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); + + 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++; + + 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; + + 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]+1-((int)calagedate %12)/12.)] += weight[i]; + freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i]; + } + } + } + } + + for(i=(int)agemin; i <= (int)agemax+3; 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]; + } + + 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( i <= (int) agemax){ + if(pos>=1.e-5){ + probs[i][jk][j1]= pp[jk]/pos; + } + } + } + + } + } + } + + + free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); + free_vector(pp,1,nlstate); + } /* End of Freq */ /************* Waves Concatenation ***************/ @@ -1270,11 +1376,11 @@ void concatwav(int wav[], int **dh, int /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; double sum=0., jmean=0.;*/ -int j, k=0,jk, ju, jl; - double sum=0.; -jmin=1e+5; - jmax=-1; -jmean=0.; + int j, k=0,jk, ju, jl; + double sum=0.; + jmin=1e+5; + jmax=-1; + jmean=0.; for(i=1; i<=imx; i++){ mi=0; m=firstpass; @@ -1309,9 +1415,9 @@ jmean=0.; if(j==0) j=1; /* Survives at least one month after exam */ k=k+1; if (j >= jmax) jmax=j; - else if (j <= jmin)jmin=j; + if (j <= jmin) jmin=j; sum=sum+j; - if (j<0) printf("j=%d num=%d ",j,i); + /* if (j<10) printf("j=%d num=%d ",j,i); */ } } else{ @@ -1319,6 +1425,7 @@ jmean=0.; k=k+1; if (j >= jmax) jmax=j; else if (j <= jmin)jmin=j; + /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ sum=sum+j; } jk= j/stepm; @@ -1335,7 +1442,7 @@ jmean=0.; } jmean=sum/k; printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean); -} + } /*********** Tricode ****************************/ void tricode(int *Tvar, int **nbcode, int imx) { @@ -1373,7 +1480,7 @@ void tricode(int *Tvar, int **nbcode, in for (k=0; k<19; k++) Ndum[k]=0; - for (i=1; i<=ncovmodel; i++) { + for (i=1; i<=ncovmodel-2; i++) { ij=Tvar[i]; Ndum[ij]++; } @@ -1448,7 +1555,7 @@ void varevsij(char fileres[], double *** double **dnewm,**doldm; int i, j, nhstepm, hstepm, h; int k, cptcode; - double *xp; + double *xp; double **gp, **gm; double ***gradg, ***trgradg; double ***p3mat; @@ -1484,6 +1591,12 @@ void varevsij(char fileres[], double *** } 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]; + } + for(j=1; j<= nlstate; j++){ for(h=0; h<=nhstepm; h++){ for(i=1, gp[h][j]=0.;i<=nlstate;i++) @@ -1495,12 +1608,19 @@ 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]; + } + for(j=1; j<= nlstate; j++){ for(h=0; h<=nhstepm; h++){ for(i=1, gm[h][j]=0.;i<=nlstate;i++) gm[h][j] += prlim[i][i]*p3mat[i][j][h]; } } + for(j=1; j<= nlstate; j++) for(h=0; h<=nhstepm; h++){ gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; @@ -1630,7 +1750,110 @@ void varprevlim(char fileres[], double * } +/************ Variance of one-step probabilities ******************/ +void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij) +{ + int i, j; + int k=0, cptcode; + double **dnewm,**doldm; + double *xp; + double *gp, *gm; + double **gradg, **trgradg; + double age,agelim, cov[NCOVMAX]; + int theta; + char fileresprob[FILENAMELENGTH]; + + strcpy(fileresprob,"prob"); + strcat(fileresprob,fileres); + if((ficresprob=fopen(fileresprob,"w"))==NULL) { + printf("Problem with resultfile: %s\n", fileresprob); + } + printf("Computing variance of one-step probabilities: result on file '%s' \n",fileresprob); + + xp=vector(1,npar); + dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); + doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath)); + + cov[1]=1; + for (age=bage; age<=fage; age ++){ + cov[2]=age; + gradg=matrix(1,npar,1,9); + trgradg=matrix(1,9,1,npar); + gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath)); + gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath)); + + for(theta=1; theta <=npar; theta++){ + for(i=1; i<=npar; i++) + xp[i] = x[i] + (i==theta ?delti[theta]:0); + + pmij(pmmij,cov,ncovmodel,xp,nlstate); + + k=0; + for(i=1; i<= (nlstate+ndeath); i++){ + for(j=1; j<=(nlstate+ndeath);j++){ + k=k+1; + gp[k]=pmmij[i][j]; + } + } + + for(i=1; i<=npar; i++) + xp[i] = x[i] - (i==theta ?delti[theta]:0); + + + pmij(pmmij,cov,ncovmodel,xp,nlstate); + k=0; + for(i=1; i<=(nlstate+ndeath); i++){ + for(j=1; j<=(nlstate+ndeath);j++){ + k=k+1; + gm[k]=pmmij[i][j]; + } + } + + for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++) + gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta]; + } + + for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++) + for(theta=1; theta <=npar; theta++) + trgradg[j][theta]=gradg[theta][j]; + + matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov); + matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg); + + pmij(pmmij,cov,ncovmodel,x,nlstate); + + k=0; + for(i=1; i<=(nlstate+ndeath); i++){ + for(j=1; j<=(nlstate+ndeath);j++){ + k=k+1; + gm[k]=pmmij[i][j]; + } + } + + /*printf("\n%d ",(int)age); + for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){ + + + printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); + }*/ + + fprintf(ficresprob,"\n%d ",(int)age); + + for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){ + if (i== 2) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]); +if (i== 4) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]); + } + + free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath)); + free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath)); + free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); + free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); +} + free_vector(xp,1,npar); +fclose(ficresprob); + exit(0); +} /***********************************************/ /**************** Main Program *****************/ @@ -1653,18 +1876,22 @@ int main() char line[MAXLINE], linepar[MAXLINE]; char title[MAXLINE]; char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH]; - char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH]; + char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH]; char filerest[FILENAMELENGTH]; char fileregp[FILENAMELENGTH]; + char popfile[FILENAMELENGTH]; char path[80],pathc[80],pathcd[80],pathtot[80],model[20]; int firstobs=1, lastobs=10; int sdeb, sfin; /* Status at beginning and end */ int c, h , cpt,l; int ju,jl, mi; int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij; - int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; - + 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; double ftolpl=FTOL; double **prlim; @@ -1677,20 +1904,27 @@ int main() double ***eij, ***vareij; 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; + char version[80]="Imach version 64b, May 2001, INED-EUROREVES "; char *alph[]={"a","a","b","c","d","e"}, str[4]; + char z[1]="c", occ; #include #include char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80]; + /* long total_usecs; struct timeval start_time, end_time; gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ - printf("\nIMACH, Version 0.64b"); + printf("\nIMACH, Version 0.7"); printf("\nEnter the parameter file name: "); #ifdef windows @@ -1738,7 +1972,15 @@ split(pathtot, path,optionfile); fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt,model); fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt,model); - +while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + puts(line); + fputs(line,ficparo); + } + ungetc(c,ficpar); + + covar=matrix(0,NCOVMAX,1,n); cptcovn=0; if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; @@ -1772,7 +2014,8 @@ split(pathtot, path,optionfile); fprintf(ficparo,"\n"); } - npar= (nlstate+ndeath-1)*nlstate*ncovmodel; + npar= (nlstate+ndeath-1)*nlstate*ncovmodel; + p=param[1][1]; /* Reads comments: lines beginning with '#' */ @@ -1862,7 +2105,7 @@ split(pathtot, path,optionfile); tab=ivector(1,NCOVMAX); ncodemax=ivector(1,8); - i=1; + i=1; while (fgets(line, MAXLINE, fic) != NULL) { if ((i >= firstobs) && (i <=lastobs)) { @@ -1884,16 +2127,26 @@ split(pathtot, path,optionfile); cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra); } num[i]=atol(stra); - - /*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]));*/ + + /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){ + 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])); ij=ij+1;}*/ i=i+1; } } - - /*scanf("%d",i);*/ + /* printf("ii=%d", ij); + scanf("%d",i);*/ imx=i-1; /* Number of individuals */ + /* for (i=1; i<=imx; i++){ + 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]));*/ + /* Calculation of the number of parameter from char model*/ Tvar=ivector(1,15); Tprod=ivector(1,15); @@ -1975,6 +2228,13 @@ split(pathtot, path,optionfile); } /*-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++) + if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){ + anint[m][i]=9999; + s[m][i]=-1; + } for (i=1; i<=imx; i++) { agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); @@ -2031,8 +2291,8 @@ printf("Total number of individuals= %d, free_imatrix(outcome,1,maxwav+1,1,n); free_vector(moisnais,1,n); free_vector(annais,1,n); - free_matrix(mint,1,maxwav,1,n); - free_matrix(anint,1,maxwav,1,n); + /* free_matrix(mint,1,maxwav,1,n); + free_matrix(anint,1,maxwav,1,n);*/ free_vector(moisdc,1,n); free_vector(andc,1,n); @@ -2064,26 +2324,18 @@ printf("Total number of individuals= %d, } } } - - - /*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'. */ - freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax); + + pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */ - + /* For Powell, parameters are in a vector p[] starting at p[1] so we point p on param[1][1] so that p[1] maps on param[1][1][1] */ p=param[1][1]; /* *(*(*(param +1)+1)+0) */ @@ -2093,8 +2345,9 @@ printf("Total number of individuals= %d, } /*--------- results files --------------*/ - fprintf(ficres,"\ntitle=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, mle,weightopt,model); - + fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, mle,weightopt,model); + + jk=1; fprintf(ficres,"# Parameters\n"); printf("# Parameters\n"); @@ -2135,7 +2388,7 @@ printf("Total number of individuals= %d, fprintf(ficres,"\n"); } } - } + } k=1; fprintf(ficres,"# Covariance\n"); @@ -2171,11 +2424,53 @@ printf("Total number of individuals= %d, fage = agemax; } - fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); + fprintf(ficres,"# agemin agemax for life expectancy.\n"); + fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage); + fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage); + + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + puts(line); + fputs(line,ficparo); + } + ungetc(c,ficpar); + + fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mob_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav); + fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mob_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); + fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mob_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); + + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + puts(line); + fputs(line,ficparo); + } + ungetc(c,ficpar); + - -/*------------ gnuplot -------------*/ + dateprev1=anprev1+mprev1/12.+jprev1/365.; + dateprev2=anprev2+mprev2/12.+jprev2/365.; + + fscanf(ficpar,"pop_based=%d\n",&popbased); + fprintf(ficparo,"pop_based=%d\n",popbased); + fprintf(ficres,"pop_based=%d\n",popbased); + + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + puts(line); + fputs(line,ficparo); + } + ungetc(c,ficpar); + fscanf(ficpar,"popforecast=%d popfile=%s starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2); +fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2); +fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2); + + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2); + + /*------------ gnuplot -------------*/ chdir(pathcd); if((ficgp=fopen("graph.plt","w"))==NULL) { printf("Problem with file graph.gp");goto end; @@ -2279,7 +2574,7 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\" fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1); fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1); } - } + } /* proba elementaires */ for(i=1,jk=1; i <=nlstate; i++){ @@ -2338,19 +2633,15 @@ ij=1; } fclose(ficgp); + /* end gnuplot */ chdir(path); - free_matrix(agev,1,maxwav,1,imx); + free_ivector(wav,1,imx); free_imatrix(dh,1,lastpass-firstpass+1,1,imx); - free_imatrix(mw,1,lastpass-firstpass+1,1,imx); - - free_imatrix(s,1,maxwav+1,1,n); - - + free_imatrix(mw,1,lastpass-firstpass+1,1,imx); free_ivector(num,1,n); free_vector(agedc,1,n); - free_vector(weight,1,n); /*free_matrix(covar,1,NCOVMAX,1,n);*/ fclose(ficparo); fclose(ficres); @@ -2381,7 +2672,7 @@ chdir(path); printf("Problem with %s \n",optionfilehtm);goto end; } - fprintf(fichtm,"
    Imach, Version 0.64b
    + fprintf(fichtm,"
      Imach, Version 0.7
      Titre=%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
      @@ -2395,7 +2686,9 @@ Interval (in months) between two waves: - 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

      ",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres); + - Standard deviation of stationary prevalences: vpl%s
      + - Prevalences and population forecasting: f%s
      +
      ",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres); fprintf(fichtm,"
    • Graphs

      • "); @@ -2480,6 +2773,7 @@ fclose(fichtm); } } fclose(ficrespl); + /*------------- h Pij x at various ages ------------*/ strcpy(filerespij,"pij"); strcat(filerespij,fileres); @@ -2489,7 +2783,7 @@ fclose(fichtm); printf("Computing pij: result on file '%s' \n", filerespij); stepsize=(int) (stepm+YEARM-1)/YEARM; - if (stepm<=24) stepsize=2; + /*if (stepm<=24) stepsize=2;*/ agelim=AGESUP; hstepm=stepsize*YEARM; /* Every year of age */ @@ -2528,8 +2822,157 @@ fclose(fichtm); } } + /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/ + 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); + for (agedeb=bage+3; agedeb<=fage-2; agedeb++) + for (i=1; i<=nlstate;i++) + for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++) + mobaverage[(int)agedeb][i][cptcod]=0.; + + for (agedeb=bage+4; agedeb<=fage; agedeb++){ + for (i=1; i<=nlstate;i++){ + for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ + for (cpt=0;cpt<=4;cpt++){ + mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod]; + } + mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5; + } + } + } + } + + 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=(bage-((int)calagedate %12)/12.); agedeb--){ /* If stepm=6 months */ + nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); + nhstepm = nhstepm/hstepm; + /*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,stepm,hstepm,nhstepm);*/ + + 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]; + /* fprintf(ficresf," p3=%.3f p=%.3f ", 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 */ + 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"); @@ -2589,6 +3032,11 @@ fclose(fichtm); epj=vector(1,nlstate+1); for(age=bage; age <=fage ;age++){ prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); + if (popbased==1) { + for(i=1; i<=nlstate;i++) + prlim[i][i]=probs[(int)age][i][k]; + } + fprintf(ficrest," %.0f",age); for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ for(i=1, epj[j]=0.;i <=nlstate;i++) { @@ -2608,6 +3056,9 @@ fclose(fichtm); } } + + + fclose(ficreseij); fclose(ficresvij); fclose(ficrest); @@ -2653,25 +3104,28 @@ strcpy(fileresvpl,"vpl"); free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); - + free_matrix(matcov,1,npar,1,npar); free_vector(delti,1,npar); free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); - printf("End of Imach\n"); + if(erreur >0) + printf("End of Imach with error %d\n",erreur); + else printf("End of Imach\n"); /* gettimeofday(&end_time, (struct timezone*)0);*/ /* after time */ /* printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);*/ /*printf("Total time was %d uSec.\n", total_usecs);*/ /*------ End -----------*/ + end: #ifdef windows chdir(pathcd); #endif - /*system("wgnuplot graph.plt");*/ - system("../gp37mgw/wgnuplot graph.plt"); + + system("..\\gp37mgw\\wgnuplot graph.plt"); #ifdef windows while (z[0] != 'q') {