--- imach096d/src/imach.c 2002/02/20 17:02:08 1.13 +++ imach096d/src/imach.c 2002/02/20 17:15:02 1.17 @@ -73,6 +73,7 @@ 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 */ @@ -84,7 +85,7 @@ double jmean; /* Mean space between 2 wa 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,*ficresf; -FILE *ficgp, *fichtm,*ficresprob; +FILE *ficgp, *fichtm,*ficresprob,*ficpop; FILE *ficreseij; char filerese[FILENAMELENGTH]; FILE *ficresvij; @@ -1149,10 +1150,10 @@ 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, int fprev1,int lprev1,double **mint,double **anint) { /* Some frequencies */ - int i, m, jk, k1, i1, j1, bool, z1,z2,j; + int i, m, jk, k1, k2,i1, j1, bool, z1,z2,j; double ***freq; /* Frequencies */ double *pp; double pos; @@ -1160,7 +1161,7 @@ void freqsummary(char fileres[], int ag char fileresp[FILENAMELENGTH]; pp=vector(1,nlstate); - probs= ma3x(1,130 ,1,8, 1,8); + probs= ma3x(1,130 ,1,8, 1,8); strcpy(fileresp,"p"); strcat(fileresp,fileres); if((ficresp=fopen(fileresp,"w"))==NULL) { @@ -1191,13 +1192,16 @@ void freqsummary(char fileres[], int ag bool=0; } if (bool==1) { - for(m=firstpass; m<=lastpass-1; m++){ + for(m=fprev1; m<=lprev1; m++){ + k2=anint[m][i]+(mint[m][i]/12.); + if ((k2>=1984) && (k2<=1988.5)) { 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 (cptcovn>0) { fprintf(ficresp, "\n#********** Variable "); @@ -1215,7 +1219,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++) @@ -1225,10 +1229,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++){ @@ -1262,6 +1268,83 @@ void freqsummary(char fileres[], int ag } /* End of Freq */ +/************ Prevalence ********************/ +void prevalence(int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax, int fprev1,int lprev1) +{ /* Some frequencies */ + + int i, m, jk, k1, i1, j1, bool, z1,z2,j; + double ***freq; /* Frequencies */ + double *pp; + double pos; + + pp=vector(1,nlstate); + probs= ma3x(1,130 ,1,8, 1,8); + + 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=fprev1; m<=lprev1; 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]; + } + } + } + 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 ***************/ void concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int firstpass, int lastpass, int imx, int nlstate, int stepm) @@ -1492,6 +1575,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++) @@ -1503,12 +1592,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]; @@ -1767,15 +1863,18 @@ int main() 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, fprev, lprev ,fprevfore=1, lprevfore=1,nforecast,popforecast=0; int hstepm, nhstepm; + int *popage; + double bage, fage, age, agelim, agebase; double ftolpl=FTOL; double **prlim; @@ -1788,7 +1887,8 @@ int main() double ***eij, ***vareij; double **varpl; /* Variances of prevalence limits by age */ double *epj, vepp; - double kk1; + double kk1, kk2; + double *popeffectif,*popcount; char version[80]="Imach version 64b, May 2001, INED-EUROREVES "; char *alph[]={"a","a","b","c","d","e"}, str[4]; @@ -1804,7 +1904,7 @@ int main() 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 @@ -1852,7 +1952,39 @@ 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); + + fscanf(ficpar,"fprevalence=%d lprevalence=%d pop_based=%d\n",&fprev,&lprev,&popbased); + fprintf(ficparo,"fprevalence=%d lprevalence=%d pop_based=%d\n",fprev,lprev,popbased); + + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + puts(line); + fputs(line,ficparo); + } + ungetc(c,ficpar); + + fscanf(ficpar,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",&fprevfore,&lprevfore,&nforecast,&mobilav); + fprintf(ficparo,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",fprevfore,lprevfore,nforecast,mobilav); + + +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\n",&popforecast,popfile); + covar=matrix(0,NCOVMAX,1,n); cptcovn=0; if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; @@ -2014,8 +2146,8 @@ split(pathtot, path,optionfile); 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++) 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++) 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); @@ -2161,8 +2293,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); @@ -2194,21 +2326,15 @@ 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 */ + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprev, lprev,mint,anint); + + free_matrix(mint,1,maxwav,1,n); + free_matrix(anint,1,maxwav,1,n); + + 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 */ @@ -2223,8 +2349,10 @@ 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); + fprintf(ficres,"fprevalence=%d lprevalence=%d pop_based=%d\n",fprev,lprev,popbased); + fprintf(ficres,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",fprevfore,lprevfore,nforecast,mobilav); + jk=1; fprintf(ficres,"# Parameters\n"); printf("# Parameters\n"); @@ -2470,17 +2598,12 @@ ij=1; fclose(ficgp); 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); @@ -2511,7 +2634,7 @@ chdir(path); printf("Problem with %s \n",optionfilehtm);goto end; } - fprintf(fichtm,"