From d9b59135e8250410b140c6a0aa87c068a9be6e2b Mon Sep 17 00:00:00 2001 From: "N. Brouard" Date: Sun, 16 May 2004 15:05:56 +0000 Subject: [PATCH] New version 0.97 . First attempt to estimate force of mortality directly from the data i.e. without the need of knowing the health state at each age, but using a Gompertz model: log u =a + b*age . This is the basic analysis of mortality and should be done before any other analysis, in order to test if the mortality estimated from the cross-longitudinal survey is different from the mortality estimated from other sources like vital statistic data. MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit The same imach parameter file can be used but the option for mle should be -3. Agn�s, who wrote this part of the code, tried to keep most of the former routines in order to include the new code within the former code. The output is very simple: only an estimate of the intercept and of the slope with 95% confident intervals. Current limitations: A) Even if you enter covariates, i.e. with the model= V1+V2 equation for example, the programm does only estimate a unique global model without covariates. B) There is no computation of Life Expectancy nor Life Table. --- src/imach.c | 1462 +++++++++++++++++++++++++++++---------------------- 1 file changed, 822 insertions(+), 640 deletions(-) diff --git a/src/imach.c b/src/imach.c index 5f33f63..303c80b 100644 --- a/src/imach.c +++ b/src/imach.c @@ -1,6 +1,10 @@ /* $Id$ $State$ $Log$ + Revision 1.97 2004/02/20 13:25:42 lievre + Version 0.96d. Population forecasting command line is (temporarily) + suppressed. + Revision 1.96 2003/07/15 15:38:55 brouard * imach.c (Repository): Errors in subdirf, 2, 3 while printing tmpout is rewritten within the same printf. Workaround: many printfs. @@ -197,6 +201,7 @@ #define YEARM 12. /* Number of months per year */ #define AGESUP 130 #define AGEBASE 40 +#define AGEGOMP 10. /* Minimal age for Gompertz adjustment */ #ifdef unix #define DIRSEPARATOR '/' #define ODIRSEPARATOR '\\' @@ -208,7 +213,7 @@ /* $Id$ */ /* $State$ */ -char version[]="Imach version 0.96d, February 2004, INED-EUROREVES "; +char version[]="Imach version 0.70, May 2004, INED-EUROREVES "; char fullversion[]="$Revision$ $Date$"; int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ int nvar; @@ -238,6 +243,7 @@ int globpr; /* Global variable for printing or not */ double fretone; /* Only one call to likelihood */ long ipmx; /* Number of contributions */ double sw; /* Sum of weights */ +char filerespow[FILENAMELENGTH]; char fileresilk[FILENAMELENGTH]; /* File of individual contributions to the likelihood */ FILE *ficresilk; FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor; @@ -300,9 +306,10 @@ static double maxarg1,maxarg2; static double sqrarg; #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 :sqrarg*sqrarg) #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;} +int agegomp= AGEGOMP; int imx; -int stepm; +int stepm=1; /* Stepm, step in month: minimum step interpolation*/ int estepm; @@ -310,9 +317,10 @@ int estepm; int m,nb; long *num; -int firstpass=0, lastpass=4,*cod, *ncodemax, *Tage; +int firstpass=0, lastpass=4,*cod, *ncodemax, *Tage,*cens; double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint; double **pmmij, ***probs; +double *ageexmed,*agecens; double dateintmean=0; double *weight; @@ -825,9 +833,10 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret last_time=curr_time; (void) gettimeofday(&curr_time,&tzp); printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout); - fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); + /* fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); - for (i=1;i<=n;i++) { + */ + for (i=1;i<=n;i++) { printf(" %d %.12f",i, p[i]); fprintf(ficlog," %d %.12lf",i, p[i]); fprintf(ficrespow," %.12lf", p[i]); @@ -1512,7 +1521,7 @@ void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, doub double **xi; double fret; double fretone; /* Only one call to likelihood */ - char filerespow[FILENAMELENGTH]; + /* char filerespow[FILENAMELENGTH];*/ xi=matrix(1,npar,1,npar); for (i=1;i<=npar;i++) for (j=1;j<=npar;j++) @@ -1547,11 +1556,11 @@ void hesscov(double **matcov, double p[], int npar, double delti[], double ftolh int i, j,jk; int *indx; - double hessii(double p[], double delta, int theta, double delti[]); - double hessij(double p[], double delti[], int i, int j); + double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar); + double hessij(double p[], double delti[], int i, int j,double (*func)(double []),int npar); void lubksb(double **a, int npar, int *indx, double b[]) ; void ludcmp(double **a, int npar, int *indx, double *d) ; - + double gompertz(double p[]); hess=matrix(1,npar,1,npar); printf("\nCalculation of the hessian matrix. Wait...\n"); @@ -1559,9 +1568,11 @@ void hesscov(double **matcov, double p[], int npar, double delti[], double ftolh for (i=1;i<=npar;i++){ printf("%d",i);fflush(stdout); fprintf(ficlog,"%d",i);fflush(ficlog); - hess[i][i]=hessii(p,ftolhess,i,delti); - /*printf(" %f ",p[i]);*/ - /*printf(" %lf ",hess[i][i]);*/ + + hess[i][i]=hessii(p,ftolhess,i,delti,func,npar); + + /* printf(" %f ",p[i]); + printf(" %lf %lf %lf",hess[i][i],ftolhess,delti[i]);*/ } for (i=1;i<=npar;i++) { @@ -1569,7 +1580,8 @@ void hesscov(double **matcov, double p[], int npar, double delti[], double ftolh if (j>i) { printf(".%d%d",i,j);fflush(stdout); fprintf(ficlog,".%d%d",i,j);fflush(ficlog); - hess[i][j]=hessij(p,delti,i,j); + hess[i][j]=hessij(p,delti,i,j,func,npar); + hess[j][i]=hess[i][j]; /*printf(" %lf ",hess[i][j]);*/ } @@ -1640,14 +1652,14 @@ void hesscov(double **matcov, double p[], int npar, double delti[], double ftolh } /*************** hessian matrix ****************/ -double hessii( double x[], double delta, int theta, double delti[]) +double hessii(double x[], double delta, int theta, double delti[], double (*func)(double []), int npar) { int i; int l=1, lmax=20; double k1,k2; double p2[NPARMAX+1]; double res; - double delt, delts, nkhi=10.,nkhif=1., khi=1.e-4; + double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; double fx; int k=0,kmax=10; double l1; @@ -1687,7 +1699,7 @@ double hessii( double x[], double delta, int theta, double delti[]) } -double hessij( double x[], double delti[], int thetai,int thetaj) +double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) { int i; int l=1, l1, lmax=20; @@ -3863,6 +3875,82 @@ void prwizard(int ncovmodel, int nlstate, int ndeath, char model[], FILE *ficpa } /* end itimes */ } /* end of prwizard */ +/******************* Gompertz Likelihood ******************************/ +double gompertz(double x[]) +{ + double A,B,L=0.0,sump=0.,num=0.; + int i,n=0; /* n is the size of the sample */ + for (i=0;i<=imx-1 ; i++) { + sump=sump+weight[i]; + sump=sump+1; + num=num+1; + } + + + /* for (i=1; i<=imx; i++) + if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/ + + for (i=0;i<=imx-1 ; i++) + { + if (cens[i]==1 & wav[i]>1) + A=-x[1]/(x[2])* + (exp(x[2]/YEARM*(agecens[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12))); + + if (cens[i]==0 & wav[i]>1) + A=-x[1]/(x[2])* + (exp(x[2]/YEARM*(agedc[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12))) + +log(x[1]/YEARM)+x[2]/YEARM*(agedc[i]*12-agegomp*12)+log(YEARM); + + if (wav[i]>1 & agecens[i]>15) { + L=L+A*weight[i]; + /* printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/ + } + } + + /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/ + + return -2*L*num/sump; +} + +/******************* Printing html file ***********/ +void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \ + int lastpass, int stepm, int weightopt, char model[],\ + int imx, double p[],double **matcov){ + int i; + + fprintf(fichtm,""); + fflush(fichtm); +} + +/******************* Gnuplot file **************/ +void printinggnuplotmort(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ + + char dirfileres[132],optfileres[132]; + int m,cpt,k1,i,k,j,jk,k2,k3,ij,l; + int ng; + + + /*#ifdef windows */ + fprintf(ficgp,"cd \"%s\" \n",pathc); + /*#endif */ + + + strcpy(dirfileres,optionfilefiname); + strcpy(optfileres,"vpl"); + fprintf(ficgp,"set out \"graphmort.png\"\n "); + fprintf(ficgp,"set xlabel \"Age\"\n set ylabel \"Force of mortality (per year)\" \n "); + fprintf(ficgp, "set ter png small\n set log y\n"); + fprintf(ficgp, "set size 0.65,0.65\n"); + fprintf(ficgp,"plot [%d:100] %lf*exp(%lf*(x-%d))",agegomp,p[1],p[2],agegomp); + +} + + /***********************************************/ @@ -3876,6 +3964,7 @@ int main(int argc, char *argv[]) int jj, ll, li, lj, lk, imk; int numlinepar=0; /* Current linenumber of parameter file */ int itimes; + int NDIM=2; char ca[32], cb[32], cc[32]; /* FILE *fichtm; *//* Html File */ @@ -3919,9 +4008,9 @@ int main(int argc, char *argv[]) double *epj, vepp; double kk1, kk2; double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000; - + double **ximort; char *alph[]={"a","a","b","c","d","e"}, str[4]; - + int *dcwave; char z[1]="c", occ; @@ -4077,163 +4166,174 @@ int main(int argc, char *argv[]) ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ - + npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ + + delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); + delti=delti3[1][1]; + /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ if(mle==-1){ /* Print a wizard for help writing covariance matrix */ prwizard(ncovmodel, nlstate, ndeath, model, ficparo); printf(" You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso); fprintf(ficlog," You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso); + free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); fclose (ficparo); fclose (ficlog); exit(0); } - /* Read guess parameters */ - /* Reads comments: lines beginning with '#' */ - while((c=getc(ficpar))=='#' && c!= EOF){ - ungetc(c,ficpar); - fgets(line, MAXLINE, ficpar); - numlinepar++; - puts(line); - fputs(line,ficparo); - fputs(line,ficlog); + else if(mle==-3) { + prwizard(ncovmodel, nlstate, ndeath, model, ficparo); + printf(" You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso); + fprintf(ficlog," You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso); + param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); + matcov=matrix(1,npar,1,npar); } - ungetc(c,ficpar); - - param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); - for(i=1; i <=nlstate; i++){ - j=0; - for(jj=1; jj <=nlstate+ndeath; jj++){ - if(jj==i) continue; - j++; - fscanf(ficpar,"%1d%1d",&i1,&j1); - if ((i1 != i) && (j1 != j)){ - printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); - exit(1); - } - fprintf(ficparo,"%1d%1d",i1,j1); - if(mle==1) - printf("%1d%1d",i,j); - fprintf(ficlog,"%1d%1d",i,j); - for(k=1; k<=ncovmodel;k++){ - fscanf(ficpar," %lf",¶m[i][j][k]); - if(mle==1){ - printf(" %lf",param[i][j][k]); - fprintf(ficlog," %lf",param[i][j][k]); - } - else - fprintf(ficlog," %lf",param[i][j][k]); - fprintf(ficparo," %lf",param[i][j][k]); - } - fscanf(ficpar,"\n"); + else{ + /* Read guess parameters */ + /* Reads comments: lines beginning with '#' */ + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); numlinepar++; - if(mle==1) - printf("\n"); - fprintf(ficlog,"\n"); - fprintf(ficparo,"\n"); + puts(line); + fputs(line,ficparo); + fputs(line,ficlog); } - } - fflush(ficlog); + ungetc(c,ficpar); + + param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); + for(i=1; i <=nlstate; i++){ + j=0; + for(jj=1; jj <=nlstate+ndeath; jj++){ + if(jj==i) continue; + j++; + fscanf(ficpar,"%1d%1d",&i1,&j1); + if ((i1 != i) && (j1 != j)){ + printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); + exit(1); + } + fprintf(ficparo,"%1d%1d",i1,j1); + if(mle==1) + printf("%1d%1d",i,j); + fprintf(ficlog,"%1d%1d",i,j); + for(k=1; k<=ncovmodel;k++){ + fscanf(ficpar," %lf",¶m[i][j][k]); + if(mle==1){ + printf(" %lf",param[i][j][k]); + fprintf(ficlog," %lf",param[i][j][k]); + } + else + fprintf(ficlog," %lf",param[i][j][k]); + fprintf(ficparo," %lf",param[i][j][k]); + } + fscanf(ficpar,"\n"); + numlinepar++; + if(mle==1) + printf("\n"); + fprintf(ficlog,"\n"); + fprintf(ficparo,"\n"); + } + } + fflush(ficlog); - npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ - p=param[1][1]; - - /* Reads comments: lines beginning with '#' */ - while((c=getc(ficpar))=='#' && c!= EOF){ + p=param[1][1]; + + /* Reads comments: lines beginning with '#' */ + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + numlinepar++; + puts(line); + fputs(line,ficparo); + fputs(line,ficlog); + } ungetc(c,ficpar); - fgets(line, MAXLINE, ficpar); - numlinepar++; - puts(line); - fputs(line,ficparo); - fputs(line,ficlog); - } - ungetc(c,ficpar); - delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); - /* delti=vector(1,npar); *//* Scale of each paramater (output from hesscov) */ - for(i=1; i <=nlstate; i++){ - for(j=1; j <=nlstate+ndeath-1; j++){ - fscanf(ficpar,"%1d%1d",&i1,&j1); - if ((i1-i)*(j1-j)!=0){ - printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); - exit(1); - } - printf("%1d%1d",i,j); - fprintf(ficparo,"%1d%1d",i1,j1); - fprintf(ficlog,"%1d%1d",i1,j1); - for(k=1; k<=ncovmodel;k++){ - fscanf(ficpar,"%le",&delti3[i][j][k]); - printf(" %le",delti3[i][j][k]); - fprintf(ficparo," %le",delti3[i][j][k]); - fprintf(ficlog," %le",delti3[i][j][k]); + for(i=1; i <=nlstate; i++){ + for(j=1; j <=nlstate+ndeath-1; j++){ + fscanf(ficpar,"%1d%1d",&i1,&j1); + if ((i1-i)*(j1-j)!=0){ + printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); + exit(1); + } + printf("%1d%1d",i,j); + fprintf(ficparo,"%1d%1d",i1,j1); + fprintf(ficlog,"%1d%1d",i1,j1); + for(k=1; k<=ncovmodel;k++){ + fscanf(ficpar,"%le",&delti3[i][j][k]); + printf(" %le",delti3[i][j][k]); + fprintf(ficparo," %le",delti3[i][j][k]); + fprintf(ficlog," %le",delti3[i][j][k]); + } + fscanf(ficpar,"\n"); + numlinepar++; + printf("\n"); + fprintf(ficparo,"\n"); + fprintf(ficlog,"\n"); } - fscanf(ficpar,"\n"); - numlinepar++; - printf("\n"); - fprintf(ficparo,"\n"); - fprintf(ficlog,"\n"); } - } - fflush(ficlog); + fflush(ficlog); - delti=delti3[1][1]; + delti=delti3[1][1]; - /* free_ma3x(delti3,1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */ /* Hasn't to to freed here otherwise delti is no more allocated */ + /* free_ma3x(delti3,1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */ /* Hasn't to to freed here otherwise delti is no more allocated */ - /* Reads comments: lines beginning with '#' */ - while((c=getc(ficpar))=='#' && c!= EOF){ + /* Reads comments: lines beginning with '#' */ + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + numlinepar++; + puts(line); + fputs(line,ficparo); + fputs(line,ficlog); + } ungetc(c,ficpar); - fgets(line, MAXLINE, ficpar); - numlinepar++; - puts(line); - fputs(line,ficparo); - fputs(line,ficlog); - } - ungetc(c,ficpar); - matcov=matrix(1,npar,1,npar); - for(i=1; i <=npar; i++){ - fscanf(ficpar,"%s",&str); - if(mle==1) - printf("%s",str); - fprintf(ficlog,"%s",str); - fprintf(ficparo,"%s",str); - for(j=1; j <=i; j++){ - fscanf(ficpar," %le",&matcov[i][j]); - if(mle==1){ - printf(" %.5le",matcov[i][j]); + matcov=matrix(1,npar,1,npar); + for(i=1; i <=npar; i++){ + fscanf(ficpar,"%s",&str); + if(mle==1) + printf("%s",str); + fprintf(ficlog,"%s",str); + fprintf(ficparo,"%s",str); + for(j=1; j <=i; j++){ + fscanf(ficpar," %le",&matcov[i][j]); + if(mle==1){ + printf(" %.5le",matcov[i][j]); + } + fprintf(ficlog," %.5le",matcov[i][j]); + fprintf(ficparo," %.5le",matcov[i][j]); } - fprintf(ficlog," %.5le",matcov[i][j]); - fprintf(ficparo," %.5le",matcov[i][j]); + fscanf(ficpar,"\n"); + numlinepar++; + if(mle==1) + printf("\n"); + fprintf(ficlog,"\n"); + fprintf(ficparo,"\n"); } - fscanf(ficpar,"\n"); - numlinepar++; + for(i=1; i <=npar; i++) + for(j=i+1;j<=npar;j++) + matcov[i][j]=matcov[j][i]; + if(mle==1) printf("\n"); fprintf(ficlog,"\n"); - fprintf(ficparo,"\n"); - } - for(i=1; i <=npar; i++) - for(j=i+1;j<=npar;j++) - matcov[i][j]=matcov[j][i]; - - if(mle==1) - printf("\n"); - fprintf(ficlog,"\n"); - - fflush(ficlog); - - /*-------- Rewriting paramater file ----------*/ - strcpy(rfileres,"r"); /* "Rparameterfile */ - strcat(rfileres,optionfilefiname); /* Parameter file first name*/ - strcat(rfileres,"."); /* */ - strcat(rfileres,optionfilext); /* Other files have txt extension */ - if((ficres =fopen(rfileres,"w"))==NULL) { - printf("Problem writing new parameter file: %s\n", fileres);goto end; - fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end; - } - fprintf(ficres,"#%s\n",version); + fflush(ficlog); + + /*-------- Rewriting parameter file ----------*/ + strcpy(rfileres,"r"); /* "Rparameterfile */ + strcat(rfileres,optionfilefiname); /* Parameter file first name*/ + strcat(rfileres,"."); /* */ + strcat(rfileres,optionfilext); /* Other files have txt extension */ + if((ficres =fopen(rfileres,"w"))==NULL) { + printf("Problem writing new parameter file: %s\n", fileres);goto end; + fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end; + } + fprintf(ficres,"#%s\n",version); + } /* End of mle != -3 */ + /*-------- data file ----------*/ if((fic=fopen(datafile,"r"))==NULL) { printf("Problem with datafile: %s\n", datafile);goto end; @@ -4488,6 +4588,7 @@ int main(int argc, char *argv[]) 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); + agegomp=(int)agemin; free_vector(severity,1,maxwav); free_imatrix(outcome,1,maxwav+1,1,n); free_vector(moisnais,1,n); @@ -4541,7 +4642,10 @@ int main(int argc, char *argv[]) /*------------ gnuplot -------------*/ strcpy(optionfilegnuplot,optionfilefiname); + if(mle==-3) + strcat(optionfilegnuplot,"-mort"); strcat(optionfilegnuplot,".gp"); + if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { printf("Problem with file %s",optionfilegnuplot); } @@ -4554,6 +4658,8 @@ int main(int argc, char *argv[]) /*--------- index.htm --------*/ strcpy(optionfilehtm,optionfilefiname); /* Main html file */ + if(mle==-3) + strcat(optionfilehtm,"-mort"); strcat(optionfilehtm,".htm"); if((fichtm=fopen(optionfilehtm,"w"))==NULL) { printf("Problem with %s \n",optionfilehtm), exit(0); @@ -4611,567 +4717,643 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ p=param[1][1]; /* *(*(*(param +1)+1)+0) */ globpr=0; /* To get the number ipmx of contributions and the sum of weights*/ - likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ - printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); - for (k=1; k<=npar;k++) - printf(" %d %8.5f",k,p[k]); - printf("\n"); - globpr=1; /* to print the contributions */ - likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ - printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); - for (k=1; k<=npar;k++) - printf(" %d %8.5f",k,p[k]); - printf("\n"); - if(mle>=1){ /* Could be 1 or 2 */ - mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); + if (mle==-3){ + ximort=matrix(1,NDIM,1,NDIM); + cens=ivector(1,n); + ageexmed=vector(1,n); + agecens=vector(1,n); + dcwave=ivector(1,n); + + for (i=1; i<=imx; i++){ + dcwave[i]=-1; + for (j=1; j<=lastpass; j++) + if (s[j][i]>nlstate) { + dcwave[i]=j; + /* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/ + break; + } + } + + for (i=1; i<=imx; i++) { + if (wav[i]>0){ + ageexmed[i]=agev[mw[1][i]][i]; + j=wav[i];agecens[i]=1.; + if (ageexmed[i]>1 & wav[i]>0) agecens[i]=agev[mw[j][i]][i]; + cens[i]=1; + + if (ageexmed[i]<1) cens[i]=-1; + if (agedc[i]< AGESUP & agedc[i]>1 & dcwave[i]>firstpass & dcwave[i]<=lastpass) cens[i]=0 ; + } + else cens[i]=-1; + } + + for (i=1;i<=NDIM;i++) { + for (j=1;j<=NDIM;j++) + ximort[i][j]=(i == j ? 1.0 : 0.0); + } + + p[1]=0.1; p[2]=0.1; + /*printf("%lf %lf", p[1], p[2]);*/ + + + printf("Powell\n"); fprintf(ficlog,"Powell\n"); + strcpy(filerespow,"pow-mort"); + strcat(filerespow,fileres); + if((ficrespow=fopen(filerespow,"w"))==NULL) { + printf("Problem with resultfile: %s\n", filerespow); + fprintf(ficlog,"Problem with resultfile: %s\n", filerespow); } + fprintf(ficrespow,"# Powell\n# iter -2*LL"); + /* for (i=1;i<=nlstate;i++) + for(j=1;j<=nlstate+ndeath;j++) + if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); + */ + fprintf(ficrespow,"\n"); + + powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); + fclose(ficrespow); - /*--------- results files --------------*/ - fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model); - + hesscov(matcov, p, NDIM,delti, 1e-4, gompertz); - fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); - printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); - fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); - for(i=1,jk=1; i <=nlstate; i++){ - for(k=1; k <=(nlstate+ndeath); k++){ - if (k != i) { - printf("%d%d ",i,k); - fprintf(ficlog,"%d%d ",i,k); - fprintf(ficres,"%1d%1d ",i,k); - for(j=1; j <=ncovmodel; j++){ - printf("%f ",p[jk]); - fprintf(ficlog,"%f ",p[jk]); - fprintf(ficres,"%f ",p[jk]); - jk++; - } - printf("\n"); - fprintf(ficlog,"\n"); - fprintf(ficres,"\n"); + for(i=1; i <=NDIM; i++) + for(j=i+1;j<=NDIM;j++) + matcov[i][j]=matcov[j][i]; + + printf("\nCovariance matrix\n "); + for(i=1; i <=NDIM; i++) { + for(j=1;j<=NDIM;j++){ + printf("%f ",matcov[i][j]); } + printf("\n "); } - } - if(mle!=0){ - /* Computing hessian and covariance matrix */ - ftolhess=ftol; /* Usually correct */ - hesscov(matcov, p, npar, delti, ftolhess, func); - } - fprintf(ficres,"# Scales (for hessian or gradient estimation)\n"); - printf("# Scales (for hessian or gradient estimation)\n"); - fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n"); - for(i=1,jk=1; i <=nlstate; i++){ - for(j=1; j <=nlstate+ndeath; j++){ - if (j!=i) { - fprintf(ficres,"%1d%1d",i,j); - printf("%1d%1d",i,j); - fprintf(ficlog,"%1d%1d",i,j); - for(k=1; k<=ncovmodel;k++){ - printf(" %.5e",delti[jk]); - fprintf(ficlog," %.5e",delti[jk]); - fprintf(ficres," %.5e",delti[jk]); - jk++; + + printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp); + for (i=1;i<=NDIM;i++) + printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); + replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */ + printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); + + printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \ + stepm, weightopt,\ + model,imx,p,matcov); + } /* Endof if mle==-3 */ + + else{ /* For mle >=1 */ + + likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ + printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); + for (k=1; k<=npar;k++) + printf(" %d %8.5f",k,p[k]); + printf("\n"); + globpr=1; /* to print the contributions */ + likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ + printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); + for (k=1; k<=npar;k++) + printf(" %d %8.5f",k,p[k]); + printf("\n"); + if(mle>=1){ /* Could be 1 or 2 */ + mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); + } + + /*--------- results files --------------*/ + fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model); + + + fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); + printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); + fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); + for(i=1,jk=1; i <=nlstate; i++){ + for(k=1; k <=(nlstate+ndeath); k++){ + if (k != i) { + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + fprintf(ficres,"%1d%1d ",i,k); + for(j=1; j <=ncovmodel; j++){ + printf("%f ",p[jk]); + fprintf(ficlog,"%f ",p[jk]); + fprintf(ficres,"%f ",p[jk]); + jk++; + } + printf("\n"); + fprintf(ficlog,"\n"); + fprintf(ficres,"\n"); } - printf("\n"); - fprintf(ficlog,"\n"); - fprintf(ficres,"\n"); } } - } - - fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); - if(mle>=1) - printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); - fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); -/* # 121 Var(a12)\n\ */ -/* # 122 Cov(b12,a12) Var(b12)\n\ */ -/* # 131 Cov(a13,a12) Cov(a13,b12, Var(a13)\n\ */ -/* # 132 Cov(b13,a12) Cov(b13,b12, Cov(b13,a13) Var(b13)\n\ */ -/* # 212 Cov(a21,a12) Cov(a21,b12, Cov(a21,a13) Cov(a21,b13) Var(a21)\n\ */ -/* # 212 Cov(b21,a12) Cov(b21,b12, Cov(b21,a13) Cov(b21,b13) Cov(b21,a21) Var(b21)\n\ */ -/* # 232 Cov(a23,a12) Cov(a23,b12, Cov(a23,a13) Cov(a23,b13) Cov(a23,a21) Cov(a23,b21) Var(a23)\n\ */ -/* # 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n" */ - - -/* Just to have a covariance matrix which will be more understandable - even is we still don't want to manage dictionary of variables -*/ - for(itimes=1;itimes<=2;itimes++){ - jj=0; - for(i=1; i <=nlstate; i++){ + if(mle!=0){ + /* Computing hessian and covariance matrix */ + ftolhess=ftol; /* Usually correct */ + hesscov(matcov, p, npar, delti, ftolhess, func); + } + fprintf(ficres,"# Scales (for hessian or gradient estimation)\n"); + printf("# Scales (for hessian or gradient estimation)\n"); + fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n"); + for(i=1,jk=1; i <=nlstate; i++){ for(j=1; j <=nlstate+ndeath; j++){ - if(j==i) continue; - for(k=1; k<=ncovmodel;k++){ - jj++; - ca[0]= k+'a'-1;ca[1]='\0'; - if(itimes==1){ - if(mle>=1) - printf("#%1d%1d%d",i,j,k); - fprintf(ficlog,"#%1d%1d%d",i,j,k); - fprintf(ficres,"#%1d%1d%d",i,j,k); - }else{ - if(mle>=1) - printf("%1d%1d%d",i,j,k); - fprintf(ficlog,"%1d%1d%d",i,j,k); - fprintf(ficres,"%1d%1d%d",i,j,k); + if (j!=i) { + fprintf(ficres,"%1d%1d",i,j); + printf("%1d%1d",i,j); + fprintf(ficlog,"%1d%1d",i,j); + for(k=1; k<=ncovmodel;k++){ + printf(" %.5e",delti[jk]); + fprintf(ficlog," %.5e",delti[jk]); + fprintf(ficres," %.5e",delti[jk]); + jk++; } - ll=0; - for(li=1;li <=nlstate; li++){ - for(lj=1;lj <=nlstate+ndeath; lj++){ - if(lj==li) continue; - for(lk=1;lk<=ncovmodel;lk++){ - ll++; - if(ll<=jj){ - cb[0]= lk +'a'-1;cb[1]='\0'; - if(ll=1) - printf(" Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj); - fprintf(ficlog," Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj); - fprintf(ficres," Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj); - }else{ - if(mle>=1) - printf(" %.5e",matcov[jj][ll]); - fprintf(ficlog," %.5e",matcov[jj][ll]); - fprintf(ficres," %.5e",matcov[jj][ll]); - } - }else{ - if(itimes==1){ - if(mle>=1) - printf(" Var(%s%1d%1d)",ca,i,j); - fprintf(ficlog," Var(%s%1d%1d)",ca,i,j); - fprintf(ficres," Var(%s%1d%1d)",ca,i,j); + printf("\n"); + fprintf(ficlog,"\n"); + fprintf(ficres,"\n"); + } + } + } + + fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); + if(mle>=1) + printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); + fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); + /* # 121 Var(a12)\n\ */ + /* # 122 Cov(b12,a12) Var(b12)\n\ */ + /* # 131 Cov(a13,a12) Cov(a13,b12, Var(a13)\n\ */ + /* # 132 Cov(b13,a12) Cov(b13,b12, Cov(b13,a13) Var(b13)\n\ */ + /* # 212 Cov(a21,a12) Cov(a21,b12, Cov(a21,a13) Cov(a21,b13) Var(a21)\n\ */ + /* # 212 Cov(b21,a12) Cov(b21,b12, Cov(b21,a13) Cov(b21,b13) Cov(b21,a21) Var(b21)\n\ */ + /* # 232 Cov(a23,a12) Cov(a23,b12, Cov(a23,a13) Cov(a23,b13) Cov(a23,a21) Cov(a23,b21) Var(a23)\n\ */ + /* # 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n" */ + + + /* Just to have a covariance matrix which will be more understandable + even is we still don't want to manage dictionary of variables + */ + for(itimes=1;itimes<=2;itimes++){ + jj=0; + for(i=1; i <=nlstate; i++){ + for(j=1; j <=nlstate+ndeath; j++){ + if(j==i) continue; + for(k=1; k<=ncovmodel;k++){ + jj++; + ca[0]= k+'a'-1;ca[1]='\0'; + if(itimes==1){ + if(mle>=1) + printf("#%1d%1d%d",i,j,k); + fprintf(ficlog,"#%1d%1d%d",i,j,k); + fprintf(ficres,"#%1d%1d%d",i,j,k); + }else{ + if(mle>=1) + printf("%1d%1d%d",i,j,k); + fprintf(ficlog,"%1d%1d%d",i,j,k); + fprintf(ficres,"%1d%1d%d",i,j,k); + } + ll=0; + for(li=1;li <=nlstate; li++){ + for(lj=1;lj <=nlstate+ndeath; lj++){ + if(lj==li) continue; + for(lk=1;lk<=ncovmodel;lk++){ + ll++; + if(ll<=jj){ + cb[0]= lk +'a'-1;cb[1]='\0'; + if(ll=1) + printf(" Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj); + fprintf(ficlog," Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj); + fprintf(ficres," Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj); + }else{ + if(mle>=1) + printf(" %.5e",matcov[jj][ll]); + fprintf(ficlog," %.5e",matcov[jj][ll]); + fprintf(ficres," %.5e",matcov[jj][ll]); + } }else{ - if(mle>=1) - printf(" %.5e",matcov[jj][ll]); - fprintf(ficlog," %.5e",matcov[jj][ll]); - fprintf(ficres," %.5e",matcov[jj][ll]); + if(itimes==1){ + if(mle>=1) + printf(" Var(%s%1d%1d)",ca,i,j); + fprintf(ficlog," Var(%s%1d%1d)",ca,i,j); + fprintf(ficres," Var(%s%1d%1d)",ca,i,j); + }else{ + if(mle>=1) + printf(" %.5e",matcov[jj][ll]); + fprintf(ficlog," %.5e",matcov[jj][ll]); + fprintf(ficres," %.5e",matcov[jj][ll]); + } } } - } - } /* end lk */ - } /* end lj */ - } /* end li */ - if(mle>=1) - printf("\n"); - fprintf(ficlog,"\n"); - fprintf(ficres,"\n"); - numlinepar++; - } /* end k*/ - } /*end j */ - } /* end i */ - } /* end itimes */ - - fflush(ficlog); - fflush(ficres); - - while((c=getc(ficpar))=='#' && c!= EOF){ - ungetc(c,ficpar); - fgets(line, MAXLINE, ficpar); - puts(line); - fputs(line,ficparo); - } - ungetc(c,ficpar); - - estepm=0; - fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); - if (estepm==0 || estepm < stepm) estepm=stepm; - if (fage <= 2) { - bage = ageminpar; - fage = agemaxpar; - } - - 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 estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); - fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); - - while((c=getc(ficpar))=='#' && c!= EOF){ + } /* end lk */ + } /* end lj */ + } /* end li */ + if(mle>=1) + printf("\n"); + fprintf(ficlog,"\n"); + fprintf(ficres,"\n"); + numlinepar++; + } /* end k*/ + } /*end j */ + } /* end i */ + } /* end itimes */ + + fflush(ficlog); + fflush(ficres); + + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + puts(line); + fputs(line,ficparo); + } 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 mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav); - fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); - fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); - printf("begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); - fprintf(ficlog,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); - - while((c=getc(ficpar))=='#' && c!= EOF){ + + estepm=0; + fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); + if (estepm==0 || estepm < stepm) estepm=stepm; + if (fage <= 2) { + bage = ageminpar; + fage = agemaxpar; + } + + 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 estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); + fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); + + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + puts(line); + fputs(line,ficparo); + } ungetc(c,ficpar); - fgets(line, MAXLINE, ficpar); - puts(line); - fputs(line,ficparo); - } - ungetc(c,ficpar); - - - dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.; - dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/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){ + + fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav); + fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); + fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); + printf("begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); + fprintf(ficlog,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_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); - fgets(line, MAXLINE, ficpar); - puts(line); - fputs(line,ficparo); - } - ungetc(c,ficpar); - - fscanf(ficpar,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj); - fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); - printf("prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); - fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); - fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); - /* day and month of proj2 are not used but only year anproj2.*/ - - while((c=getc(ficpar))=='#' && c!= EOF){ + + + dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.; + dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/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); - fgets(line, MAXLINE, ficpar); - puts(line); - fputs(line,ficparo); - } - ungetc(c,ficpar); - - /* fscanf(ficpar,"popforecast=%d popfile=%s popfiledate=%lf/%lf/%lf last-popfiledate=%lf/%lf/%lf\n",&popforecast,popfile,&jpyram,&mpyram,&anpyram,&jpyram1,&mpyram1,&anpyram1); - 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,Tvaraff,nbcode, ncodemax,mint,anint);*/ - /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ - - replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */ - printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); - - printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\ - model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\ - jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); - - /*------------ free_vector -------------*/ - /* chdir(path); */ + + fscanf(ficpar,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj); + fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); + printf("prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); + fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); + fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); + /* day and month of proj2 are not used but only year anproj2.*/ + + + + /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);*/ + /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ + + replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */ + printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); + + printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\ + model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\ + jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); + + /*------------ free_vector -------------*/ + /* chdir(path); */ - free_ivector(wav,1,imx); - free_imatrix(dh,1,lastpass-firstpass+1,1,imx); - free_imatrix(bh,1,lastpass-firstpass+1,1,imx); - free_imatrix(mw,1,lastpass-firstpass+1,1,imx); - free_lvector(num,1,n); - free_vector(agedc,1,n); - /*free_matrix(covar,0,NCOVMAX,1,n);*/ - /*free_matrix(covar,1,NCOVMAX,1,n);*/ - fclose(ficparo); - fclose(ficres); - - - /*--------------- Prevalence limit (stable prevalence) --------------*/ + free_ivector(wav,1,imx); + free_imatrix(dh,1,lastpass-firstpass+1,1,imx); + free_imatrix(bh,1,lastpass-firstpass+1,1,imx); + free_imatrix(mw,1,lastpass-firstpass+1,1,imx); + free_lvector(num,1,n); + free_vector(agedc,1,n); + /*free_matrix(covar,0,NCOVMAX,1,n);*/ + /*free_matrix(covar,1,NCOVMAX,1,n);*/ + fclose(ficparo); + fclose(ficres); + + + /*--------------- Prevalence limit (stable prevalence) --------------*/ - strcpy(filerespl,"pl"); - strcat(filerespl,fileres); - if((ficrespl=fopen(filerespl,"w"))==NULL) { - printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end; - fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end; - } - printf("Computing stable prevalence: result on file '%s' \n", filerespl); - fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl); - fprintf(ficrespl,"#Stable prevalence \n"); - fprintf(ficrespl,"#Age "); - for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); - fprintf(ficrespl,"\n"); + strcpy(filerespl,"pl"); + strcat(filerespl,fileres); + if((ficrespl=fopen(filerespl,"w"))==NULL) { + printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end; + fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end; + } + printf("Computing stable prevalence: result on file '%s' \n", filerespl); + fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl); + fprintf(ficrespl,"#Stable prevalence \n"); + fprintf(ficrespl,"#Age "); + for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); + fprintf(ficrespl,"\n"); - prlim=matrix(1,nlstate,1,nlstate); - - agebase=ageminpar; - agelim=agemaxpar; - ftolpl=1.e-10; - i1=cptcoveff; - if (cptcovn < 1){i1=1;} - - for(cptcov=1,k=0;cptcov<=i1;cptcov++){ - for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ - k=k+1; - /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/ - fprintf(ficrespl,"\n#******"); - printf("\n#******"); - fprintf(ficlog,"\n#******"); - for(j=1;j<=cptcoveff;j++) { - fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); - printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); - fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); - } - fprintf(ficrespl,"******\n"); - printf("******\n"); - fprintf(ficlog,"******\n"); + prlim=matrix(1,nlstate,1,nlstate); + + agebase=ageminpar; + agelim=agemaxpar; + ftolpl=1.e-10; + i1=cptcoveff; + if (cptcovn < 1){i1=1;} + + for(cptcov=1,k=0;cptcov<=i1;cptcov++){ + for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ + k=k+1; + /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/ + fprintf(ficrespl,"\n#******"); + printf("\n#******"); + fprintf(ficlog,"\n#******"); + for(j=1;j<=cptcoveff;j++) { + fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + } + fprintf(ficrespl,"******\n"); + printf("******\n"); + fprintf(ficlog,"******\n"); - for (age=agebase; age<=agelim; age++){ - prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); - fprintf(ficrespl,"%.0f ",age ); - for(j=1;j<=cptcoveff;j++) - fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); - for(i=1; i<=nlstate;i++) - fprintf(ficrespl," %.5f", prlim[i][i]); - fprintf(ficrespl,"\n"); + for (age=agebase; age<=agelim; age++){ + prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); + fprintf(ficrespl,"%.0f ",age ); + for(j=1;j<=cptcoveff;j++) + fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + for(i=1; i<=nlstate;i++) + fprintf(ficrespl," %.5f", prlim[i][i]); + fprintf(ficrespl,"\n"); + } } } - } - fclose(ficrespl); + fclose(ficrespl); - /*------------- h Pij x at various ages ------------*/ + /*------------- h Pij x at various ages ------------*/ - strcpy(filerespij,"pij"); strcat(filerespij,fileres); - if((ficrespij=fopen(filerespij,"w"))==NULL) { - printf("Problem with Pij resultfile: %s\n", filerespij);goto end; - fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij);goto end; - } - printf("Computing pij: result on file '%s' \n", filerespij); - fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij); + strcpy(filerespij,"pij"); strcat(filerespij,fileres); + if((ficrespij=fopen(filerespij,"w"))==NULL) { + printf("Problem with Pij resultfile: %s\n", filerespij);goto end; + fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij);goto end; + } + printf("Computing pij: result on file '%s' \n", filerespij); + fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij); - stepsize=(int) (stepm+YEARM-1)/YEARM; - /*if (stepm<=24) stepsize=2;*/ - - agelim=AGESUP; - hstepm=stepsize*YEARM; /* Every year of age */ - hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ - - /* hstepm=1; aff par mois*/ - - fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x "); - for(cptcov=1,k=0;cptcov<=i1;cptcov++){ - for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ - k=k+1; - fprintf(ficrespij,"\n#****** "); - for(j=1;j<=cptcoveff;j++) - fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); - fprintf(ficrespij,"******\n"); + stepsize=(int) (stepm+YEARM-1)/YEARM; + /*if (stepm<=24) stepsize=2;*/ + + agelim=AGESUP; + hstepm=stepsize*YEARM; /* Every year of age */ + hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ + + /* hstepm=1; aff par mois*/ + + fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x "); + for(cptcov=1,k=0;cptcov<=i1;cptcov++){ + for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ + k=k+1; + fprintf(ficrespij,"\n#****** "); + for(j=1;j<=cptcoveff;j++) + fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficrespij,"******\n"); - for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */ - nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ - nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ + for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */ + nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ + nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ - /* nhstepm=nhstepm*YEARM; aff par mois*/ + /* nhstepm=nhstepm*YEARM; aff par mois*/ - 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); - fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j="); - for(i=1; i<=nlstate;i++) - for(j=1; j<=nlstate+ndeath;j++) - fprintf(ficrespij," %1d-%1d",i,j); - fprintf(ficrespij,"\n"); - for (h=0; h<=nhstepm; h++){ - fprintf(ficrespij,"%d %3.f %3.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm ); + 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); + fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j="); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate+ndeath;j++) - fprintf(ficrespij," %.5f", p3mat[i][j][h]); + fprintf(ficrespij," %1d-%1d",i,j); + fprintf(ficrespij,"\n"); + for (h=0; h<=nhstepm; h++){ + fprintf(ficrespij,"%d %3.f %3.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm ); + for(i=1; i<=nlstate;i++) + for(j=1; j<=nlstate+ndeath;j++) + fprintf(ficrespij," %.5f", p3mat[i][j][h]); + fprintf(ficrespij,"\n"); + } + free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); fprintf(ficrespij,"\n"); } - free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - fprintf(ficrespij,"\n"); } } - } - varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax); + varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax); - fclose(ficrespij); + fclose(ficrespij); - probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); - for(i=1;i<=AGESUP;i++) - for(j=1;j<=NCOVMAX;j++) - for(k=1;k<=NCOVMAX;k++) - probs[i][j][k]=0.; + probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); + for(i=1;i<=AGESUP;i++) + for(j=1;j<=NCOVMAX;j++) + for(k=1;k<=NCOVMAX;k++) + probs[i][j][k]=0.; - /*---------- Forecasting ------------------*/ - /*if((stepm == 1) && (strcmp(model,".")==0)){*/ - if(prevfcast==1){ - /* if(stepm ==1){*/ + /*---------- Forecasting ------------------*/ + /*if((stepm == 1) && (strcmp(model,".")==0)){*/ + if(prevfcast==1){ + /* if(stepm ==1){*/ prevforecast(fileres, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff); /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/ -/* } */ -/* else{ */ -/* erreur=108; */ -/* printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */ -/* fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */ -/* } */ - } + /* } */ + /* else{ */ + /* erreur=108; */ + /* printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */ + /* fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */ + /* } */ + } - /*---------- Health expectancies and variances ------------*/ + /*---------- Health expectancies and variances ------------*/ - strcpy(filerest,"t"); - strcat(filerest,fileres); - if((ficrest=fopen(filerest,"w"))==NULL) { - printf("Problem with total LE resultfile: %s\n", filerest);goto end; - fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end; - } - printf("Computing Total LEs with variances: file '%s' \n", filerest); - fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest); + strcpy(filerest,"t"); + strcat(filerest,fileres); + if((ficrest=fopen(filerest,"w"))==NULL) { + printf("Problem with total LE resultfile: %s\n", filerest);goto end; + fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end; + } + printf("Computing Total LEs with variances: file '%s' \n", filerest); + fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest); - strcpy(filerese,"e"); - strcat(filerese,fileres); - if((ficreseij=fopen(filerese,"w"))==NULL) { - printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0); - fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0); - } - printf("Computing Health Expectancies: result on file '%s' \n", filerese); - fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese); - - strcpy(fileresv,"v"); - strcat(fileresv,fileres); - if((ficresvij=fopen(fileresv,"w"))==NULL) { - printf("Problem with variance resultfile: %s\n", fileresv);exit(0); - fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0); - } - printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); - fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); + strcpy(filerese,"e"); + strcat(filerese,fileres); + if((ficreseij=fopen(filerese,"w"))==NULL) { + printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0); + fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0); + } + printf("Computing Health Expectancies: result on file '%s' \n", filerese); + fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese); + + strcpy(fileresv,"v"); + strcat(fileresv,fileres); + if((ficresvij=fopen(fileresv,"w"))==NULL) { + printf("Problem with variance resultfile: %s\n", fileresv);exit(0); + fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0); + } + printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); + fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); - /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */ - prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); - /* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\ -ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass); - */ + /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */ + prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); + /* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\ + ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass); + */ - if (mobilav!=0) { - mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ - fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); - printf(" Error in movingaverage mobilav=%d\n",mobilav); + if (mobilav!=0) { + mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ + fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); + printf(" Error in movingaverage mobilav=%d\n",mobilav); + } } - } - for(cptcov=1,k=0;cptcov<=i1;cptcov++){ - for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ - k=k+1; - fprintf(ficrest,"\n#****** "); - for(j=1;j<=cptcoveff;j++) - fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); - fprintf(ficrest,"******\n"); - - fprintf(ficreseij,"\n#****** "); - for(j=1;j<=cptcoveff;j++) - fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); - fprintf(ficreseij,"******\n"); - - fprintf(ficresvij,"\n#****** "); - for(j=1;j<=cptcoveff;j++) - fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); - fprintf(ficresvij,"******\n"); - - eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); - oldm=oldms;savm=savms; - evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov); + for(cptcov=1,k=0;cptcov<=i1;cptcov++){ + for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ + k=k+1; + fprintf(ficrest,"\n#****** "); + for(j=1;j<=cptcoveff;j++) + fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficrest,"******\n"); + + fprintf(ficreseij,"\n#****** "); + for(j=1;j<=cptcoveff;j++) + fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficreseij,"******\n"); + + fprintf(ficresvij,"\n#****** "); + for(j=1;j<=cptcoveff;j++) + fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficresvij,"******\n"); + + eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); + oldm=oldms;savm=savms; + evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov); - vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); - oldm=oldms;savm=savms; - varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav); - if(popbased==1){ - varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav); - } + vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); + oldm=oldms;savm=savms; + varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav); + if(popbased==1){ + varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav); + } - fprintf(ficrest,"#Total LEs with variances: e.. (std) "); - for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); - fprintf(ficrest,"\n"); - - epj=vector(1,nlstate+1); - for(age=bage; age <=fage ;age++){ - prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); - if (popbased==1) { - if(mobilav ==0){ - for(i=1; i<=nlstate;i++) - prlim[i][i]=probs[(int)age][i][k]; - }else{ /* mobilav */ - for(i=1; i<=nlstate;i++) - prlim[i][i]=mobaverage[(int)age][i][k]; + fprintf(ficrest,"#Total LEs with variances: e.. (std) "); + for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); + fprintf(ficrest,"\n"); + + epj=vector(1,nlstate+1); + for(age=bage; age <=fage ;age++){ + prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); + if (popbased==1) { + if(mobilav ==0){ + for(i=1; i<=nlstate;i++) + prlim[i][i]=probs[(int)age][i][k]; + }else{ /* mobilav */ + for(i=1; i<=nlstate;i++) + prlim[i][i]=mobaverage[(int)age][i][k]; + } } - } - fprintf(ficrest," %4.0f",age); - for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ - for(i=1, epj[j]=0.;i <=nlstate;i++) { - epj[j] += prlim[i][i]*eij[i][j][(int)age]; - /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ + fprintf(ficrest," %4.0f",age); + for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ + for(i=1, epj[j]=0.;i <=nlstate;i++) { + epj[j] += prlim[i][i]*eij[i][j][(int)age]; + /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ + } + epj[nlstate+1] +=epj[j]; } - epj[nlstate+1] +=epj[j]; - } - for(i=1, vepp=0.;i <=nlstate;i++) - for(j=1;j <=nlstate;j++) - vepp += vareij[i][j][(int)age]; - fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); - for(j=1;j <=nlstate;j++){ - fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); + for(i=1, vepp=0.;i <=nlstate;i++) + for(j=1;j <=nlstate;j++) + vepp += vareij[i][j][(int)age]; + fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); + for(j=1;j <=nlstate;j++){ + fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); + } + fprintf(ficrest,"\n"); } - fprintf(ficrest,"\n"); + free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); + free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); + free_vector(epj,1,nlstate+1); } - free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); - free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); - free_vector(epj,1,nlstate+1); } - } - free_vector(weight,1,n); - free_imatrix(Tvard,1,15,1,2); - free_imatrix(s,1,maxwav+1,1,n); - free_matrix(anint,1,maxwav,1,n); - free_matrix(mint,1,maxwav,1,n); - free_ivector(cod,1,n); - free_ivector(tab,1,NCOVMAX); - fclose(ficreseij); - fclose(ficresvij); - fclose(ficrest); - fclose(ficpar); + free_vector(weight,1,n); + free_imatrix(Tvard,1,15,1,2); + free_imatrix(s,1,maxwav+1,1,n); + free_matrix(anint,1,maxwav,1,n); + free_matrix(mint,1,maxwav,1,n); + free_ivector(cod,1,n); + free_ivector(tab,1,NCOVMAX); + fclose(ficreseij); + fclose(ficresvij); + fclose(ficrest); + fclose(ficpar); - /*------- Variance of stable prevalence------*/ + /*------- Variance of stable prevalence------*/ - strcpy(fileresvpl,"vpl"); - strcat(fileresvpl,fileres); - if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { - printf("Problem with variance of stable prevalence resultfile: %s\n", fileresvpl); - exit(0); - } - printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl); - - for(cptcov=1,k=0;cptcov<=i1;cptcov++){ - for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ - k=k+1; - fprintf(ficresvpl,"\n#****** "); - for(j=1;j<=cptcoveff;j++) - fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); - fprintf(ficresvpl,"******\n"); + strcpy(fileresvpl,"vpl"); + strcat(fileresvpl,fileres); + if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { + printf("Problem with variance of stable prevalence resultfile: %s\n", fileresvpl); + exit(0); + } + printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl); + + for(cptcov=1,k=0;cptcov<=i1;cptcov++){ + for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ + k=k+1; + fprintf(ficresvpl,"\n#****** "); + for(j=1;j<=cptcoveff;j++) + fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficresvpl,"******\n"); - varpl=matrix(1,nlstate,(int) bage, (int) fage); - oldm=oldms;savm=savms; - varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k); - free_matrix(varpl,1,nlstate,(int) bage, (int)fage); + varpl=matrix(1,nlstate,(int) bage, (int) fage); + oldm=oldms;savm=savms; + varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k); + free_matrix(varpl,1,nlstate,(int) bage, (int)fage); + } } - } - fclose(ficresvpl); + fclose(ficresvpl); - /*---------- End : free ----------------*/ - free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); - 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); + /*---------- End : free ----------------*/ + if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); + + } /* mle==-3 arrives here for freeing */ + free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); + 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(covar,0,NCOVMAX,1,n); - free_matrix(matcov,1,npar,1,npar); - /*free_vector(delti,1,npar);*/ - free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); - free_matrix(agev,1,maxwav,1,imx); - free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); - if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); - free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); - - free_ivector(ncodemax,1,8); - free_ivector(Tvar,1,15); - free_ivector(Tprod,1,15); - free_ivector(Tvaraff,1,15); - free_ivector(Tage,1,15); - free_ivector(Tcode,1,100); + free_matrix(covar,0,NCOVMAX,1,n); + free_matrix(matcov,1,npar,1,npar); + /*free_vector(delti,1,npar);*/ + free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); + free_matrix(agev,1,maxwav,1,imx); + free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); + + free_ivector(ncodemax,1,8); + free_ivector(Tvar,1,15); + free_ivector(Tprod,1,15); + free_ivector(Tvaraff,1,15); + free_ivector(Tage,1,15); + free_ivector(Tcode,1,100); + fflush(fichtm); fflush(ficgp); -- 2.43.0