--- imach/src/imach.c 2002/07/25 07:37:44 1.57 +++ imach/src/imach.c 2002/07/26 12:29:55 1.58 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.57 2002/07/25 07:37:44 lievre Exp $ +/* $Id: imach.c,v 1.58 2002/07/26 12:29:55 lievre Exp $ Interpolated Markov Chain Short summary of the programme: @@ -1284,7 +1284,7 @@ void freqsummary(char fileres[], int ag if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) bool=0; } - if (bool==1) { + if (bool==1){ for(m=firstpass; m<=lastpass; m++){ k2=anint[m][i]+(mint[m][i]/12.); if ((k2>=dateprev1) && (k2<=dateprev2)) { @@ -1575,31 +1575,35 @@ void concatwav(int wav[], int **dh, int /*********** Tricode ****************************/ void tricode(int *Tvar, int **nbcode, int imx) { - int Ndum[20],ij=1, k, j, i; + + int Ndum[20],ij=1, k, j, i, maxncov=19; int cptcode=0; cptcoveff=0; - for (k=0; k<19; k++) Ndum[k]=0; + for (k=0; k cptcode) cptcode=ij; + if (ij > cptcode) cptcode=ij; /* getting the maximum of covariable + Tvar[j]. If V=sex and male is 0 and + female is 1, then cptcode=1.*/ } for (i=0; i<=cptcode; i++) { - if(Ndum[i]!=0) ncodemax[j]++; + if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariates. In fact ncodemax[j]=2 (dichotom. variables) but it can be more */ } - ij=1; - + ij=1; for (i=1; i<=ncodemax[j]; i++) { - for (k=0; k<=19; k++) { + for (k=0; k<= maxncov; k++) { if (Ndum[k] != 0) { nbcode[Tvar[j]][ij]=k; + /* store the modality in an array. k is a modality. If we have model=V1+V1*sex then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ ij++; } @@ -1608,22 +1612,23 @@ void tricode(int *Tvar, int **nbcode, in } } - for (k=0; k<19; k++) Ndum[k]=0; + for (k=0; k< maxncov; k++) Ndum[k]=0; - for (i=1; i<=ncovmodel-2; i++) { + for (i=1; i<=ncovmodel-2; i++) { + /* Listing of all covariables in staement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ ij=Tvar[i]; - Ndum[ij]++; + Ndum[ij]++; } ij=1; - for (i=1; i<=10; i++) { + for (i=1; i<= maxncov; i++) { if((Ndum[i]!=0) && (i<=ncovcol)){ - Tvaraff[ij]=i; + Tvaraff[ij]=i; /*For printing */ ij++; } } - cptcoveff=ij-1; + cptcoveff=ij-1; /*Number of simple covariates*/ } /*********** Health Expectancies ****************/ @@ -1824,7 +1829,7 @@ void varevsij(char optionfilefiname[], d char fileresprobmorprev[FILENAMELENGTH]; if(popbased==1){ - if(mobilav==1) + if(mobilav!=0) strcpy(digitp,"-populbased-mobilav-"); else strcpy(digitp,"-populbased-nomobil-"); } @@ -2773,14 +2778,20 @@ m=pow(2,cptcoveff); int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav){ int i, cpt, cptcod; + int modcovmax =1; int mobilavrange, mob; double age; + + modcovmax=2*cptcoveff;/* Max number of modalities. We suppose + a covariate has 2 modalities */ + if (cptcovn<1) modcovmax=1; /* At least 1 pass */ + if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){ if(mobilav==1) mobilavrange=5; /* default */ else mobilavrange=mobilav; for (age=bage; age<=fage; age++) for (i=1; i<=nlstate;i++) - for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++) + for (cptcod=1;cptcod<=modcovmax;cptcod++) mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod]; /* We keep the original values on the extreme ages bage, fage and for fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2 @@ -2789,7 +2800,7 @@ int movingaverage(double ***probs, doubl for (mob=3;mob <=mobilavrange;mob=mob+2){ for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ for (i=1; i<=nlstate;i++){ - for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ + for (cptcod=1;cptcod<=modcovmax;cptcod++){ mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod]; for (cpt=1;cpt<=(mob-1)/2;cpt++){ mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod]; @@ -2817,7 +2828,7 @@ prevforecast(char fileres[], double anpr char fileresf[FILENAMELENGTH]; agelim=AGESUP; -calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; + calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); @@ -3209,10 +3220,10 @@ while((c=getc(ficpar))=='#' && c!= EOF){ covar=matrix(0,NCOVMAX,1,n); - cptcovn=0; + cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/ if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; - ncovmodel=2+cptcovn; + ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ /* Read guess parameters */ @@ -3250,7 +3261,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ fprintf(ficparo,"\n"); } - npar= (nlstate+ndeath-1)*nlstate*ncovmodel; + npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ p=param[1][1]; @@ -3409,12 +3420,12 @@ while((c=getc(ficpar))=='#' && c!= EOF){ Tvard=imatrix(1,15,1,2); Tage=ivector(1,15); - if (strlen(model) >1){ + if (strlen(model) >1){ /* If there is at least 1 covariate */ j=0, j1=0, k1=1, k2=1; - j=nbocc(model,'+'); - j1=nbocc(model,'*'); - cptcovn=j+1; - cptcovprod=j1; + j=nbocc(model,'+'); /* j=Number of '+' */ + j1=nbocc(model,'*'); /* j1=Number of '*' */ + cptcovn=j+1; + cptcovprod=j1; /*Number of products */ strcpy(modelsav,model); if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){ @@ -3423,6 +3434,8 @@ while((c=getc(ficpar))=='#' && c!= EOF){ goto end; } + /* This loop fill the array Tvar from the string 'model'.*/ + for(i=(j+1); i>=1;i--){ cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */ if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyze it */ @@ -3472,11 +3485,15 @@ while((c=getc(ficpar))=='#' && c!= EOF){ } /* end of loop + */ } /* end model */ + /*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products. + If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/ + /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]); printf("cptcovprod=%d ", cptcovprod); fprintf(ficlog,"cptcovprod=%d ", cptcovprod); - scanf("%d ",i);*/ - fclose(fic); + + scanf("%d ",i); + fclose(fic);*/ /* if(mle==1){*/ if (weightopt != 1) { /* Maximisation without weights*/ @@ -3547,8 +3564,8 @@ while((c=getc(ficpar))=='#' && c!= EOF){ } } -printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); - fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); + printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); + fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); free_vector(severity,1,maxwav); free_imatrix(outcome,1,maxwav+1,1,n); @@ -3567,13 +3584,15 @@ printf("Total number of individuals= %d, /* Concatenates waves */ concatwav(wav, dh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); + /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ Tcode=ivector(1,100); nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); ncodemax[1]=1; if (cptcovn > 0) tricode(Tvar,nbcode,imx); - codtab=imatrix(1,100,1,10); + codtab=imatrix(1,100,1,10); /* Cross tabulation to get the order of + the estimations*/ h=0; m=pow(2,cptcoveff); @@ -3770,7 +3789,8 @@ while((c=getc(ficpar))=='#' && c!= EOF){ 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,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); + /*------------ gnuplot -------------*/ strcpy(optionfilegnuplot,optionfilefiname); strcat(optionfilegnuplot,".gp"); @@ -3975,8 +3995,11 @@ Interval (in months) between two waves: } printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); + calagedate=-1; + prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); + if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){