--- imach/src/imach.c 2001/05/02 17:50:24 1.7 +++ imach/src/imach.c 2001/05/02 17:54:31 1.8 @@ -68,8 +68,7 @@ int nvar; -static int cptcov; -int cptcovn, cptcovage=0, cptcoveff=0; +int cptcovn, cptcovage=0, cptcoveff=0,cptcov; int npar=NPARMAX; int nlstate=2; /* Number of live states */ int ndeath=1; /* Number of dead states */ @@ -77,9 +76,11 @@ int ncovmodel, ncov; /* Total number int *wav; /* Number of waves for this individuual 0 is possible */ int maxwav; /* Maxim number of waves */ +int jmin, jmax; /* min, max spacing between 2 waves */ int mle, weightopt; int **mw; /* mw[mi][i] is number of the mi wave for this individual */ int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */ +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; @@ -666,7 +667,6 @@ double **prevalim(double **prlim, int nl 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]);*/ out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); @@ -841,34 +841,31 @@ double func( double *x) /* We are differentiating ll according to initial status */ /* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ /*for(i=1;i0) { for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) bool=0; + if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) + bool=0; } if (bool==1) { for(m=firstpass; m<=lastpass-1; m++){ @@ -1199,8 +1198,8 @@ void freqsummary(char fileres[], int ag 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"); @@ -1268,9 +1267,14 @@ void concatwav(int wav[], int **dh, int */ int i, mi, m; - int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; -float sum=0.; + /* 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.; for(i=1; i<=imx; i++){ mi=0; m=firstpass; @@ -1301,10 +1305,16 @@ float sum=0.; else{ if (s[mw[mi+1][i]][i] > nlstate) { j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); - if(j=0) j=1; /* Survives at least one month after exam */ + /*if ((j<0) || (j>28)) printf("j=%d num=%d ",j,i);*/ + 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; + sum=sum+j; } else{ j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); + /*if ((j<0) || (j>28)) printf("j=%d num=%d ",j,i);*/ k=k+1; if (j >= jmax) jmax=j; else if (j <= jmin)jmin=j; @@ -1322,7 +1332,8 @@ float sum=0.; } } } - printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,sum/k); + 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) @@ -1338,45 +1349,43 @@ void tricode(int *Tvar, int **nbcode, in for (i=1; i<=imx; i++) { ij=(int)(covar[Tvar[j]][i]); Ndum[ij]++; + /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/ if (ij > cptcode) cptcode=ij; } - /*printf("cptcode=%d cptcovn=%d ",cptcode,cptcovn);*/ for (i=0; i<=cptcode; i++) { if(Ndum[i]!=0) ncodemax[j]++; } ij=1; + for (i=1; i<=ncodemax[j]; i++) { for (k=0; k<=19; k++) { if (Ndum[k] != 0) { nbcode[Tvar[j]][ij]=k; - /* printf("ij=%d ",nbcode[Tvar[2]][1]);*/ ij++; } if (ij > ncodemax[j]) break; } } } - for (i=1; i<=10; i++) { + + for (k=0; k<19; k++) Ndum[k]=0; + + for (i=1; i<=ncovmodel; i++) { ij=Tvar[i]; Ndum[ij]++; } + ij=1; - for (i=1; i<=cptcovn; i++) { + for (i=1; i<=10; i++) { if((Ndum[i]!=0) && (i<=ncov)){ - Tvaraff[i]=ij; - ij++; + Tvaraff[ij]=i; + ij++; } } - for (j=1; j<=(cptcovn+2*cptcovprod); j++) { - if ((Tvar[j]>= cptcoveff) && (Tvar[j] <=ncov)) cptcoveff=Tvar[j]; - /*printf("j=%d %d\n",j,Tvar[j]);*/ - } - - /* printf("cptcoveff=%d Tvaraff=%d %d\n",cptcoveff, Tvaraff[1],Tvaraff[2]); - scanf("%d",i);*/ + cptcoveff=ij-1; } /*********** Health Expectancies ****************/ @@ -1630,7 +1639,7 @@ void varprevlim(char fileres[], double * int main() { - int i,j, k, n=MAXN,iter,m,size,cptcode, aaa, cptcod; + int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod; double agedeb, agefin,hf; double agemin=1.e20, agemax=-1.e20; @@ -1729,13 +1738,9 @@ split(pathtot, path,optionfile); 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); - covar=matrix(0,NCOVMAX,1,n); - if (strlen(model)<=1) cptcovn=0; - else { - j=0; - j=nbocc(model,'+'); - cptcovn=j+1; - } + covar=matrix(0,NCOVMAX,1,n); + cptcovn=0; + if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; ncovmodel=2+cptcovn; nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ @@ -1902,90 +1907,65 @@ split(pathtot, path,optionfile); cptcovn=j+1; cptcovprod=j1; + strcpy(modelsav,model); - if (j==0) { - if (j1==0){ - cutv(stra,strb,modelsav,'V'); - Tvar[1]=atoi(strb); - } - else if (j1==1) { - cutv(stra,strb,modelsav,'*'); - Tage[1]=1; cptcovage++; - if (strcmp(stra,"age")==0) { + if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){ + printf("Error. Non available option model=%s ",model); + goto end; + } + + for(i=(j+1); i>=1;i--){ + cutv(stra,strb,modelsav,'+'); + if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); + /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ + /*scanf("%d",i);*/ + if (strchr(strb,'*')) { + cutv(strd,strc,strb,'*'); + if (strcmp(strc,"age")==0) { cptcovprod--; - cutv(strd,strc,strb,'V'); - Tvar[1]=atoi(strc); + cutv(strb,stre,strd,'V'); + Tvar[i]=atoi(stre); + cptcovage++; + Tage[cptcovage]=i; + /*printf("stre=%s ", stre);*/ } - else if (strcmp(strb,"age")==0) { + else if (strcmp(strd,"age")==0) { cptcovprod--; - cutv(strd,strc,stra,'V'); - Tvar[1]=atoi(strc); + cutv(strb,stre,strc,'V'); + Tvar[i]=atoi(stre); + cptcovage++; + Tage[cptcovage]=i; } else { - cutv(strd,strc,strb,'V'); - cutv(stre,strd,stra,'V'); - Tvar[1]=ncov+1; + cutv(strb,stre,strc,'V'); + Tvar[i]=ncov+k1; + cutv(strb,strc,strd,'V'); + Tprod[k1]=i; + Tvard[k1][1]=atoi(strc); + Tvard[k1][2]=atoi(stre); + Tvar[cptcovn+k2]=Tvard[k1][1]; + Tvar[cptcovn+k2+1]=Tvard[k1][2]; for (k=1; k<=lastobs;k++) - covar[ncov+1][k]=covar[atoi(strc)][k]*covar[atoi(strd)][k]; + covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k]; + k1++; + k2=k2+2; } - /*printf("%s %s %s\n", stra,strb,modelsav); -printf("%d ",Tvar[1]); -scanf("%d",i);*/ - } - } - else { - for(i=j; i>=1;i--){ - cutv(stra,strb,modelsav,'+'); - /*printf("%s %s %s\n", stra,strb,modelsav); - scanf("%d",i);*/ - if (strchr(strb,'*')) { - cutv(strd,strc,strb,'*'); - if (strcmp(strc,"age")==0) { - cptcovprod--; - cutv(strb,stre,strd,'V'); - Tvar[i+1]=atoi(stre); - cptcovage++; - Tage[cptcovage]=i+1; - printf("stre=%s ", stre); - } - else if (strcmp(strd,"age")==0) { - cptcovprod--; - cutv(strb,stre,strc,'V'); - Tvar[i+1]=atoi(stre); - cptcovage++; - Tage[cptcovage]=i+1; - } - else { - cutv(strb,stre,strc,'V'); - Tvar[i+1]=ncov+k1; - cutv(strb,strc,strd,'V'); - Tprod[k1]=i+1; - Tvard[k1][1]=atoi(strc); - Tvard[k1][2]=atoi(stre); - Tvar[cptcovn+k2]=Tvard[k1][1]; - Tvar[cptcovn+k2+1]=Tvard[k1][2]; - for (k=1; k<=lastobs;k++) - covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k]; - k1++; - k2=k2+2; - } - } - else { - cutv(strd,strc,strb,'V'); - /* printf("%s %s %s", strd,strc,strb);*/ - Tvar[i+1]=atoi(strc); - } - strcpy(modelsav,stra); } - cutv(strd,strc,stra,'V'); - Tvar[1]=atoi(strc); + else { + /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ + /* scanf("%d",i);*/ + cutv(strd,strc,strb,'V'); + Tvar[i]=atoi(strc); + } + strcpy(modelsav,stra); + /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav); + scanf("%d",i);*/ } - } - /* for (i=1; i<=5; i++) - printf("i=%d %d ",i,Tvar[i]);*/ - /* printf("tvar=%d %d cptcovage=%d %d",Tvar[1],Tvar[2],cptcovage,Tage[1]);*/ - /*printf("cptcovprod=%d ", cptcovprod);*/ - /* scanf("%d ",i);*/ +} + + /*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); /* if(mle==1){*/ @@ -2003,9 +1983,11 @@ scanf("%d",i);*/ if(agedc[i]>0) if(moisdc[i]!=99 && andc[i]!=9999) agev[m][i]=agedc[i]; - else{ + else { + if (andc[i]!=9999){ printf("Warning negative age at death: %d line:%d\n",num[i],i); agev[m][i]=-1; + } } } else if(s[m][i] !=9){ /* Should no more exist */ @@ -2063,7 +2045,7 @@ printf("Total number of individuals= %d, Tcode=ivector(1,100); - nbcode=imatrix(1,nvar,1,8); + nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); ncodemax[1]=1; if (cptcovn > 0) tricode(Tvar,nbcode,imx); @@ -2396,7 +2378,12 @@ chdir(path); printf("Problem with index.htm \n");goto end; } - fprintf(fichtm,"