--- imach096d/src/imach.c 2001/05/02 17:42:45 1.5 +++ imach096d/src/imach.c 2001/05/17 16:07:14 1.11 @@ -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: @@ -54,8 +54,6 @@ #define GLOCK_ERROR_NOPATH -1 /* empty path */ #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ - - #define MAXPARM 30 /* Maximum number of parameters for the optimization */ #define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */ @@ -70,8 +68,7 @@ int nvar; -static int cptcov; -int cptcovn; +int cptcovn, cptcovage=0, cptcoveff=0,cptcov; int npar=NPARMAX; int nlstate=2; /* Number of live states */ int ndeath=1; /* Number of dead states */ @@ -79,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; @@ -93,9 +92,6 @@ FILE *ficreseij; FILE *ficresvpl; char fileresvpl[FILENAMELENGTH]; - - - #define NR_END 1 #define FREE_ARG char* #define FTOL 1.0e-10 @@ -129,19 +125,19 @@ int stepm; /* Stepm, step in month: minimum step interpolation*/ int m,nb; -int *num, firstpass=0, lastpass=4,*cod, *ncodemax; +int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage; double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint; double **pmmij; double *weight; int **s; /* Status */ double *agedc, **covar, idx; -int **nbcode, *Tcode, *Tvar, **codtab; +int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff; double ftol=FTOL; /* Tolerance for computing Max Likelihood */ double ftolhess; /* Tolerance for computing hessian */ - +/**************** split *************************/ static int split( char *path, char *dirc, char *name ) { char *s; /* pointer */ @@ -205,7 +201,7 @@ int nbocc(char *s, char occ) void cutv(char *u,char *v, char*t, char occ) { - int i,lg,j,p; + int i,lg,j,p=0; i=0; for(j=0; j<=strlen(t)-1; j++) { if((t[j]!= occ) && (t[j+1]== occ)) p=j+1; @@ -214,8 +210,8 @@ void cutv(char *u,char *v, char*t, char lg=strlen(t); for(j=0; j=(p+1))(v[j-p-1] = t[j]); @@ -652,15 +648,27 @@ double **prevalim(double **prlim, int nl for (j=1;j<=nlstate+ndeath;j++){ oldm[ii][j]=(ii==j ? 1.0 : 0.0); } - /* Even if hstepm = 1, at least one multiplication by the unit matrix */ + + cov[1]=1.; + + /* Even if hstepm = 1, at least one multiplication by the unit matrix */ for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){ newm=savm; /* Covariates have to be included here again */ - cov[1]=1.; - cov[2]=agefin; - if (cptcovn>0){ - for (k=1; k<=cptcovn;k++) {cov[2+k]=nbcode[Tvar[k]][codtab[ij][k]];/*printf("Tcode[ij]=%d nbcode=%d\n",Tcode[ij],nbcode[k][Tcode[ij]]);*/} - } + cov[2]=agefin; + + 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]);*/ + } + 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]);*/ + out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); savm=oldm; @@ -795,9 +803,13 @@ double ***hpxij(double ***po, int nhstep /* Covariates have to be included here again */ cov[1]=1.; cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM; - if (cptcovn>0){ - for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][k]]; - } + for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[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<=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("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, @@ -819,7 +831,7 @@ double ***hpxij(double ***po, int nhstep /*************** log-likelihood *************/ double func( double *x) { - int i, ii, j, k, mi, d; + int i, ii, j, k, mi, d, kk; double l, ll[NLSTATEMAX], cov[NCOVMAX]; double **out; double sw; /* Sum of weights */ @@ -829,29 +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 (k=1; k<=cptcovn;k++) cov[2+k]=covar[1+k-1][i]; - } - out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, - 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); - savm=oldm; - oldm=newm; - - + cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM; + for (kk=1; kk<=cptcovage;kk++) { + cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; + } + + out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, + 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); + savm=oldm; + oldm=newm; + + } /* end mult */ - + lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ ipmx +=1; @@ -1152,13 +1166,14 @@ void freqsummary(char fileres[], int ag freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3); j1=0; - j=cptcovn; + 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++) @@ -1167,8 +1182,9 @@ void freqsummary(char fileres[], int ag for (i=1; i<=imx; i++) { bool=1; if (cptcovn>0) { - for (z1=1; z1<=cptcovn; z1++) - if (covar[Tvar[z1]][i]!= nbcode[Tvar[z1]][codtab[j1][z1]]) bool=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-1; m++){ @@ -1180,10 +1196,10 @@ void freqsummary(char fileres[], int ag } } if (cptcovn>0) { - fprintf(ficresp, "\n#Variable"); - for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " V%d=%d",Tvar[z1],nbcode[Tvar[z1]][codtab[j1][z1]]); - } - fprintf(ficresp, "\n#"); + 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"); @@ -1251,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; @@ -1283,8 +1304,15 @@ float sum=0.; dh[mi][i]=1; else{ if (s[mw[mi+1][i]][i] > nlstate) { + if (agedc[i] < 2*AGESUP) { 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=1; /* Survives at least one month after exam */ + k=k+1; + if (j >= jmax) jmax=j; + if (j <= jmin) jmin=j; + sum=sum+j; + if (j<0) printf("j=%d num=%d ",j,i); + } } else{ j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); @@ -1305,30 +1333,35 @@ 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) { - int Ndum[80],ij, k, j, i; + int Ndum[20],ij=1, k, j, i; int cptcode=0; - for (k=0; k<79; k++) Ndum[k]=0; + cptcoveff=0; + + for (k=0; k<19; k++) Ndum[k]=0; for (k=1; k<=7; k++) ncodemax[k]=0; - - for (j=1; j<=cptcovn; j++) { + + for (j=1; j<=(cptcovn+2*cptcovprod); j++) { 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<=79; k++) { + for (k=0; k<=19; k++) { if (Ndum[k] != 0) { nbcode[Tvar[j]][ij]=k; ij++; @@ -1336,9 +1369,25 @@ void tricode(int *Tvar, int **nbcode, in if (ij > ncodemax[j]) break; } } - } + } - } + 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<=10; i++) { + if((Ndum[i]!=0) && (i<=ncov)){ + Tvaraff[ij]=i; + ij++; + } + } + + cptcoveff=ij-1; +} /*********** Health Expectancies ****************/ @@ -1591,7 +1640,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; @@ -1603,7 +1652,7 @@ int main() int *indx; char line[MAXLINE], linepar[MAXLINE]; char title[MAXLINE]; - char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH]; + char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH]; char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH]; char filerest[FILENAMELENGTH]; char fileregp[FILENAMELENGTH]; @@ -1612,7 +1661,7 @@ int main() 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; + int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij; int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; int hstepm, nhstepm; @@ -1628,7 +1677,7 @@ int main() double ***eij, ***vareij; double **varpl; /* Variances of prevalence limits by age */ double *epj, vepp; - char version[80]="Imach version 62c, May 1999, INED-EUROREVES "; + 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; @@ -1641,7 +1690,7 @@ int main() gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ - printf("\nIMACH, Version 0.64a"); + printf("\nIMACH, Version 0.64b"); printf("\nEnter the parameter file name: "); #ifdef windows @@ -1690,13 +1739,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(1,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 */ @@ -1788,7 +1833,6 @@ split(pathtot, path,optionfile); printf("\n"); - if(mle==1){ /*-------- data file ----------*/ if((ficres =fopen(fileres,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileres);goto end; @@ -1851,40 +1895,81 @@ split(pathtot, path,optionfile); imx=i-1; /* Number of individuals */ /* Calculation of the number of parameter from char model*/ - Tvar=ivector(1,8); + Tvar=ivector(1,15); + Tprod=ivector(1,15); + Tvaraff=ivector(1,15); + Tvard=imatrix(1,15,1,2); + Tage=ivector(1,15); if (strlen(model) >1){ - j=0; + j=0, j1=0, k1=1, k2=1; j=nbocc(model,'+'); + j1=nbocc(model,'*'); cptcovn=j+1; + cptcovprod=j1; + strcpy(modelsav,model); - if (j==0) { - cutv(stra,strb,modelsav,'V'); Tvar[1]=atoi(strb); + if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){ + printf("Error. Non available option model=%s ",model); + goto end; } - else { - for(i=j; i>=1;i--){ - cutv(stra,strb,modelsav,'+'); - if (strchr(strb,'*')) { - cutv(strd,strc,strb,'*'); - cutv(strb,stre,strc,'V');Tvar[i+1]=ncov+1; + + 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(strb,stre,strd,'V'); + Tvar[i]=atoi(stre); + cptcovage++; + Tage[cptcovage]=i; + /*printf("stre=%s ", stre);*/ + } + else if (strcmp(strd,"age")==0) { + cptcovprod--; + cutv(strb,stre,strc,'V'); + Tvar[i]=atoi(stre); + cptcovage++; + Tage[cptcovage]=i; + } + else { + 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(stre)][k]*covar[atoi(strc)][k]; + covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k]; + k1++; + k2=k2+2; } - else {cutv(strd,strc,strb,'V'); - 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);*/ } - } - /*printf("tvar=%d ",Tvar[1]); - 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){*/ if (weightopt != 1) { /* Maximisation without weights*/ for(i=1;i<=n;i++) weight[i]=1.0; } @@ -1899,9 +1984,11 @@ split(pathtot, path,optionfile); 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 */ @@ -1958,19 +2045,19 @@ printf("Total number of individuals= %d, concatwav(wav, dh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); -Tcode=ivector(1,100); - nbcode=imatrix(1,nvar,1,8); - ncodemax[1]=1; - if (cptcovn > 0) tricode(Tvar,nbcode,imx); - + 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); h=0; - m=pow(2,cptcovn); + m=pow(2,cptcoveff); - for(k=1;k<=cptcovn; k++){ + for(k=1;k<=cptcoveff; k++){ for(i=1; i <=(m/pow(2,k));i++){ for(j=1; j <= ncodemax[k]; j++){ - for(cpt=1; cpt <=(m/pow(2,cptcovn+1-k)); cpt++){ + for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){ h++; if (h>m) h=1;codtab[h][k]=j; } @@ -1978,13 +2065,14 @@ Tcode=ivector(1,100); } } + /*for(i=1; i <=m ;i++){ for(k=1; k <=cptcovn; k++){ - printf("i=%d k=%d %d ",i,k,codtab[i][k]); + printf("i=%d k=%d %d %d",i,k,codtab[i][k], cptcoveff); } printf("\n"); - }*/ - /*scanf("%d",i);*/ + } + scanf("%d",i);*/ /* Calculates basic frequencies. Computes observed prevalence at single age and prints on file fileres'p'. */ @@ -1999,9 +2087,10 @@ Tcode=ivector(1,100); /* 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) */ - - mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); + if(mle==1){ + mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); + } /*--------- 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); @@ -2025,10 +2114,11 @@ Tcode=ivector(1,100); } } } - + if(mle==1){ /* Computing hessian and covariance matrix */ ftolhess=ftol; /* Usually correct */ hesscov(matcov, p, npar, delti, ftolhess, func); + } fprintf(ficres,"# Scales\n"); printf("# Scales\n"); for(i=1,jk=1; i <=nlstate; i++){ @@ -2083,6 +2173,8 @@ Tcode=ivector(1,100); 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,agemax,bage,fage); + + /*------------ gnuplot -------------*/ chdir(pathcd); if((ficgp=fopen("graph.plt","w"))==NULL) { @@ -2091,7 +2183,7 @@ chdir(pathcd); #ifdef windows fprintf(ficgp,"cd \"%s\" \n",pathc); #endif -m=pow(2,cptcovn); +m=pow(2,cptcoveff); /* 1eme*/ for (cpt=1; cpt<= nlstate ; cpt ++) { @@ -2212,26 +2304,38 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\" for(k=1; k<=(nlstate+ndeath); k++) { if (k != k2){ fprintf(ficgp," exp(p%d+p%d*x",i,i+1); - - for(j=3; j <=ncovmodel; j++) - fprintf(ficgp,"+p%d*%d",k2,nbcode[Tvar[j-2]][codtab[jk][j-2]]); - fprintf(ficgp,")/(1"); - - for(k1=1; k1 <=nlstate+1; k1=k1+2){ - fprintf(ficgp,"+exp(p%d+p%d*x",k1+k3-1,k1+k3); - - for(j=3; j <=ncovmodel; j++) - fprintf(ficgp,"+p%d*%d",k2,nbcode[Tvar[j-2]][codtab[jk][j-2]]); - fprintf(ficgp,")"); +ij=1; + for(j=3; j <=ncovmodel; j++) { + if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { + fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); + ij++; + } + else + fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]); + } + fprintf(ficgp,")/(1"); + + for(k1=1; k1 <=nlstate; k1++){ + fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1); +ij=1; + for(j=3; j <=ncovmodel; j++){ + if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { + fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); + ij++; + } + else + fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]); + } + fprintf(ficgp,")"); } fprintf(ficgp,") t \"p%d%d\" ", k2,k); if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,","); - i=i+ncovmodel; + i=i+ncovmodel; } } } - fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk); - } + fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk); + } fclose(ficgp); @@ -2250,7 +2354,7 @@ chdir(path); /*free_matrix(covar,1,NCOVMAX,1,n);*/ fclose(ficparo); fclose(ficres); - } + /* }*/ /*________fin mle=1_________*/ @@ -2271,11 +2375,18 @@ chdir(path); fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage); /*--------- index.htm --------*/ - if((fichtm=fopen("index.htm","w"))==NULL) { - printf("Problem with index.htm \n");goto end; + strcpy(optionfilehtm,optionfile); + strcat(optionfilehtm,".htm"); + if((fichtm=fopen(optionfilehtm,"w"))==NULL) { + printf("Problem with %s \n",optionfilehtm);goto end; } - fprintf(fichtm,"