--- imach/src/imach.c 2003/06/24 12:34:15 1.90 +++ imach/src/imach.c 2004/09/15 17:31:30 1.102 @@ -1,6 +1,72 @@ -/* $Id: imach.c,v 1.90 2003/06/24 12:34:15 brouard Exp $ +/* $Id: imach.c,v 1.102 2004/09/15 17:31:30 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.102 2004/09/15 17:31:30 brouard + Add the possibility to read data file including tab characters. + + Revision 1.101 2004/09/15 10:38:38 brouard + Fix on curr_time + + Revision 1.100 2004/07/12 18:29:06 brouard + Add version for Mac OS X. Just define UNIX in Makefile + + Revision 1.99 2004/06/05 08:57:40 brouard + *** empty log message *** + + Revision 1.98 2004/05/16 15:05:56 brouard + 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. + + 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. + + 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. + + Revision 1.95 2003/07/08 07:54:34 brouard + * imach.c (Repository): + (Repository): Using imachwizard code to output a more meaningful covariance + matrix (cov(a12,c31) instead of numbers. + + Revision 1.94 2003/06/27 13:00:02 brouard + Just cleaning + + Revision 1.93 2003/06/25 16:33:55 brouard + (Module): On windows (cygwin) function asctime_r doesn't + exist so I changed back to asctime which exists. + (Module): Version 0.96b + + Revision 1.92 2003/06/25 16:30:45 brouard + (Module): On windows (cygwin) function asctime_r doesn't + exist so I changed back to asctime which exists. + + Revision 1.91 2003/06/25 15:30:29 brouard + * imach.c (Repository): Duplicated warning errors corrected. + (Repository): Elapsed time after each iteration is now output. It + helps to forecast when convergence will be reached. Elapsed time + is stamped in powell. We created a new html file for the graphs + concerning matrix of covariance. It has extension -cov.htm. + Revision 1.90 2003/06/24 12:34:15 brouard (Module): Some bugs corrected for windows. Also, when mle=-1 a template is output in file "or"mypar.txt with the design @@ -142,10 +208,13 @@ #include #include -#include +/* #include */ #include #include "timeval.h" +/* #include */ +/* #define _(String) gettext (String) */ + #define MAXLINE 256 #define GNUPLOTPROGRAM "gnuplot" /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ @@ -166,7 +235,8 @@ #define YEARM 12. /* Number of months per year */ #define AGESUP 130 #define AGEBASE 40 -#ifdef unix +#define AGEGOMP 10. /* Minimal age for Gompertz adjustment */ +#ifdef UNIX #define DIRSEPARATOR '/' #define ODIRSEPARATOR '\\' #else @@ -174,12 +244,12 @@ #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.90 2003/06/24 12:34:15 brouard Exp $ */ +/* $Id: imach.c,v 1.102 2004/09/15 17:31:30 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.96a, June 2003, INED-EUROREVES "; -char fullversion[]="$Revision: 1.90 $ $Date: 2003/06/24 12:34:15 $"; -int erreur; /* Error number */ +char version[]="Imach version 0.97c, September 2004, INED-EUROREVES "; +char fullversion[]="$Revision: 1.102 $ $Date: 2004/09/15 17:31:30 $"; +int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ int nvar; int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; int npar=NPARMAX; @@ -207,11 +277,12 @@ int globpr; /* Global variable for print 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; FILE *ficresprobmorprev; -FILE *fichtm; /* Html File */ +FILE *fichtm, *fichtmcov; /* Html File */ FILE *ficreseij; char filerese[FILENAMELENGTH]; FILE *ficresvij; @@ -221,18 +292,26 @@ char fileresvpl[FILENAMELENGTH]; char title[MAXLINE]; char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH]; char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH]; -char tmpout[FILENAMELENGTH]; +char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH]; char command[FILENAMELENGTH]; int outcmd=0; char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH]; -char lfileres[FILENAMELENGTH]; + char filelog[FILENAMELENGTH]; /* Log file */ char filerest[FILENAMELENGTH]; char fileregp[FILENAMELENGTH]; char popfile[FILENAMELENGTH]; -char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH]; +char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilehtmcov[FILENAMELENGTH] ; + +struct timeval start_time, end_time, curr_time, last_time, forecast_time; +struct timezone tzp; +extern int gettimeofday(); +struct tm tmg, tm, tmf, *gmtime(), *localtime(); +long time_value; +extern long time(); +char strcurr[80], strfor[80]; #define NR_END 1 #define FREE_ARG char* @@ -261,9 +340,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; @@ -271,9 +351,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; @@ -287,6 +368,9 @@ double ftolhess; /* Tolerance for comput /**************** split *************************/ static int split( char *path, char *dirc, char *name, char *ext, char *finame ) { + /* From a file name with full path (either Unix or Windows) we extract the directory (dirc) + the name of the file (name), its extension only (ext) and its first part of the name (finame) + */ char *ss; /* pointer */ int l1, l2; /* length counters */ @@ -318,12 +402,14 @@ static int split( char *path, char *dirc #endif */ ss = strrchr( name, '.' ); /* find last / */ - ss++; - strcpy(ext,ss); /* save extension */ - l1= strlen( name); - l2= strlen(ss)+1; - strncpy( finame, name, l1-l2); - finame[l1-l2]= 0; + if (ss >0){ + ss++; + strcpy(ext,ss); /* save extension */ + l1= strlen( name); + l2= strlen(ss)+1; + strncpy( finame, name, l1-l2); + finame[l1-l2]= 0; + } return( 0 ); /* we're done */ } @@ -356,8 +442,8 @@ int nbocc(char *s, char occ) void cutv(char *u,char *v, char*t, char occ) { - /* cuts string t into u and v where u is ended by char occ excluding it - and v is after occ excluding it too : ex cutv(u,v,"abcdef2ghi2j",2) + /* cuts string t into u and v where u ends before first occurence of char 'occ' + and v starts after first occurence of char 'occ' : ex cutv(u,v,"abcdef2ghi2j",'2') gives u="abcedf" and v="ghi2j" */ int i,lg,j,p=0; i=0; @@ -538,6 +624,41 @@ void free_ma3x(double ***m, long nrl, lo free((FREE_ARG)(m+nrl-NR_END)); } +/*************** function subdirf ***********/ +char *subdirf(char fileres[]) +{ + /* Caution optionfilefiname is hidden */ + strcpy(tmpout,optionfilefiname); + strcat(tmpout,"/"); /* Add to the right */ + strcat(tmpout,fileres); + return tmpout; +} + +/*************** function subdirf2 ***********/ +char *subdirf2(char fileres[], char *preop) +{ + + /* Caution optionfilefiname is hidden */ + strcpy(tmpout,optionfilefiname); + strcat(tmpout,"/"); + strcat(tmpout,preop); + strcat(tmpout,fileres); + return tmpout; +} + +/*************** function subdirf3 ***********/ +char *subdirf3(char fileres[], char *preop, char *preop2) +{ + + /* Caution optionfilefiname is hidden */ + strcpy(tmpout,optionfilefiname); + strcat(tmpout,"/"); + strcat(tmpout,preop); + strcat(tmpout,preop2); + strcat(tmpout,fileres); + return tmpout; +} + /***************** f1dim *************************/ extern int ncom; extern double *pcom,*xicom; @@ -713,6 +834,19 @@ void linmin(double p[], double xi[], int free_vector(pcom,1,n); } +char *asc_diff_time(long time_sec, char ascdiff[]) +{ + long sec_left, days, hours, minutes; + days = (time_sec) / (60*60*24); + sec_left = (time_sec) % (60*60*24); + hours = (sec_left) / (60*60) ; + sec_left = (sec_left) %(60*60); + minutes = (sec_left) /60; + sec_left = (sec_left) % (60); + sprintf(ascdiff,"%d day(s) %d hour(s) %d minute(s) %d second(s)",days, hours, minutes, sec_left); + return ascdiff; +} + /*************** powell ************************/ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret, double (*func)(double [])) @@ -723,6 +857,8 @@ void powell(double p[], double **xi, int double del,t,*pt,*ptt,*xit; double fp,fptt; double *xits; + int niterf, itmp; + pt=vector(1,n); ptt=vector(1,n); xit=vector(1,n); @@ -733,17 +869,42 @@ void powell(double p[], double **xi, int fp=(*fret); ibig=0; del=0.0; - printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret); - fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f",*iter,*fret); - fprintf(ficrespow,"%d %.12f",*iter,*fret); - for (i=1;i<=n;i++) { + 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(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); + */ + for (i=1;i<=n;i++) { printf(" %d %.12f",i, p[i]); fprintf(ficlog," %d %.12lf",i, p[i]); fprintf(ficrespow," %.12lf", p[i]); } printf("\n"); fprintf(ficlog,"\n"); - fprintf(ficrespow,"\n"); + fprintf(ficrespow,"\n");fflush(ficrespow); + if(*iter <=3){ + tm = *localtime(&curr_time.tv_sec); + strcpy(strcurr,asctime(&tm)); +/* asctime_r(&tm,strcurr); */ + forecast_time=curr_time; + itmp = strlen(strcurr); + if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */ + strcurr[itmp-1]='\0'; + printf("\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec); + fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec); + for(niterf=10;niterf<=30;niterf+=10){ + forecast_time.tv_sec=curr_time.tv_sec+(niterf-*iter)*(curr_time.tv_sec-last_time.tv_sec); + tmf = *localtime(&forecast_time.tv_sec); +/* asctime_r(&tmf,strfor); */ + strcpy(strfor,asctime(&tmf)); + itmp = strlen(strfor); + if(strfor[itmp-1]=='\n') + strfor[itmp-1]='\0'; + printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); + fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); + } + } for (i=1;i<=n;i++) { for (j=1;j<=n;j++) xit[j]=xi[j][i]; fptt=(*fret); @@ -906,57 +1067,57 @@ double **pmij(double **ps, double *cov, int i,j,j1, nc, ii, jj; for(i=1; i<= nlstate; i++){ - for(j=1; ji s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/ + for(j=1; ji s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2); */ + } + ps[i][j]=s2; } - ps[i][j]=s2; } - } /*ps[3][2]=1;*/ - - for(i=1; i<= nlstate; i++){ - s1=0; - for(j=1; j 1 the results are less biased than in previous versions. @@ -1095,15 +1256,15 @@ double func( double *x) s1=s[mw[mi][i]][i]; s2=s[mw[mi+1][i]][i]; bbh=(double)bh[mi][i]/(double)stepm; - /* bias is positive if real duration + /* bias bh is positive if real duration * is higher than the multiple of stepm and negative otherwise. */ /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ if( s2 > nlstate){ /* i.e. if s2 is a death state and if the date of death is known then the contribution to the likelihood is the probability to die between last step unit time and current - step unit time, which is also the differences between probability to die before dh - and probability to die before dh-stepm . + step unit time, which is also equal to probability to die before dh + minus probability to die before dh-stepm . In version up to 0.92 likelihood was computed as if date of death was unknown. Death was treated as any other health state: the date of the interview describes the actual state @@ -1156,31 +1317,10 @@ double func( double *x) oldm=newm; } /* end mult */ - /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ - /* But now since version 0.9 we anticipate for bias and large stepm. - * If stepm is larger than one month (smallest stepm) and if the exact delay - * (in months) between two waves is not a multiple of stepm, we rounded to - * the nearest (and in case of equal distance, to the lowest) interval but now - * we keep into memory the bias bh[mi][i] and also the previous matrix product - * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the - * probability in order to take into account the bias as a fraction of the way - * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies - * -stepm/2 to stepm/2 . - * For stepm=1 the results are the same as for previous versions of Imach. - * For stepm > 1 the results are less biased than in previous versions. - */ s1=s[mw[mi][i]][i]; s2=s[mw[mi+1][i]][i]; bbh=(double)bh[mi][i]/(double)stepm; - /* bias is positive if real duration - * is higher than the multiple of stepm and negative otherwise. - */ lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ - /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ - /*lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-+bh)*out[s1][s2])); */ /* exponential interpolation */ - /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ - /*if(lli ==000.0)*/ - /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ ipmx +=1; sw += weight[i]; ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; @@ -1207,30 +1347,10 @@ double func( double *x) oldm=newm; } /* end mult */ - /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ - /* But now since version 0.9 we anticipate for bias and large stepm. - * If stepm is larger than one month (smallest stepm) and if the exact delay - * (in months) between two waves is not a multiple of stepm, we rounded to - * the nearest (and in case of equal distance, to the lowest) interval but now - * we keep into memory the bias bh[mi][i] and also the previous matrix product - * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the - * probability in order to take into account the bias as a fraction of the way - * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies - * -stepm/2 to stepm/2 . - * For stepm=1 the results are the same as for previous versions of Imach. - * For stepm > 1 the results are less biased than in previous versions. - */ s1=s[mw[mi][i]][i]; s2=s[mw[mi+1][i]][i]; bbh=(double)bh[mi][i]/(double)stepm; - /* bias is positive if real duration - * is higher than the multiple of stepm and negative otherwise. - */ - /* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); */ /* linear interpolation */ lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ - /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ - /*if(lli ==000.0)*/ - /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ ipmx +=1; sw += weight[i]; ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; @@ -1396,35 +1516,8 @@ double funcone( double *x) return -l; } -char *subdirf(char fileres[]) -{ - - strcpy(tmpout,optionfilefiname); - strcat(tmpout,"/"); /* Add to the right */ - strcat(tmpout,fileres); - return tmpout; -} - -char *subdirf2(char fileres[], char *preop) -{ - - strcpy(tmpout,optionfilefiname); - strcat(tmpout,"/"); - strcat(tmpout,preop); - strcat(tmpout,fileres); - return tmpout; -} -char *subdirf3(char fileres[], char *preop, char *preop2) -{ - - strcpy(tmpout,optionfilefiname); - strcat(tmpout,"/"); - strcat(tmpout,preop); - strcat(tmpout,preop2); - strcat(tmpout,fileres); - return tmpout; -} +/*************** function likelione ***********/ void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*funcone)(double [])) { /* This routine should help understanding what is done with @@ -1467,7 +1560,7 @@ void mlikeli(FILE *ficres,double p[], in 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++) @@ -1502,11 +1595,11 @@ void hesscov(double **matcov, double p[] 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"); @@ -1514,9 +1607,11 @@ void hesscov(double **matcov, double p[] 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++) { @@ -1524,7 +1619,8 @@ void hesscov(double **matcov, double p[] 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]);*/ } @@ -1595,14 +1691,14 @@ void hesscov(double **matcov, double p[] } /*************** 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; @@ -1642,7 +1738,7 @@ double hessii( double x[], double delta, } -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; @@ -2036,6 +2132,7 @@ void concatwav(int wav[], int **dh, int wav[i]=mi; if(mi==0){ + nbwarn++; if(first==0){ printf("Warning! None valid information for:%ld line=%d (skipped) and may be others, see log file\n",num[i],i); first=1; @@ -2056,10 +2153,11 @@ void concatwav(int wav[], int **dh, int j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); if(j==0) j=1; /* Survives at least one month after exam */ else if(j<0){ + nberr++; printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); - j=1; /* Careful Patch */ + j=1; /* Temporary Dangerous patch */ printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n You MUST fix the contradiction between dates.\n",stepm); - printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n You MUST fix the contradiction between dates.\n",stepm); } k=k+1; @@ -2079,6 +2177,7 @@ void concatwav(int wav[], int **dh, int /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/ if(j<0){ + nberr++; printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); } @@ -2255,7 +2354,7 @@ void evsij(char fileres[], double ***eij hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ - /* Computing Variances of health expectancies */ + /* Computing Variances of health expectancies */ for(theta=1; theta <=npar; theta++){ for(i=1; i<=npar; i++){ @@ -2809,9 +2908,19 @@ void varprob(char optionfilefiname[], do fprintf(fichtm,"\n

  • Computing and drawing one step probabilities with their confidence intervals

    • \n"); fprintf(fichtm,"\n"); - fprintf(fichtm,"\n

  • Computing matrix of variance-covariance of step probabilities

    • \n"); - fprintf(fichtm,"\nWe have drawn ellipsoids of confidence around the pij, pkl to understand the covariance between two incidences. They are expressed in year-1 in order to be less dependent of stepm.
    \n"); - fprintf(fichtm,"\n
    We have drawn x'cov-1x = 4 where x is the column vector (pij,pkl). It means that if pij and pkl where uncorrelated the (2X2) matrix would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 standard deviations wide on each axis.
    When both incidences are correlated we diagonalised the inverse of the covariance matrix and made the appropriate rotation.
    \n"); + fprintf(fichtm,"\n

  • Matrix of variance-covariance of pairs of step probabilities (drawings)

    • \n",optionfilehtmcov); + fprintf(fichtmcov,"\n

    Matrix of variance-covariance of pairs of step probabilities

    \n\ + file %s
    \n",optionfilehtmcov); + fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (pij, pkl) are estimated\ +and drawn. It helps understanding how is the covariance between two incidences.\ + They are expressed in year-1 in order to be less dependent of stepm.
    \n"); + fprintf(fichtmcov,"\n
    Contour plot corresponding to x'cov-1x = 4 (where x is the column vector (pij,pkl)) are drawn. \ +It can be understood this way: if pij and pkl where uncorrelated the (2x2) matrix of covariance \ +would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 \ +standard deviations wide on each axis.
    \ + Now, if both incidences are correlated (usual case) we diagonalised the inverse of the covariance matrix\ + and made the appropriate rotation to look at the uncorrelated principal directions.
    \ +To be simple, these graphs help to understand the significativity of each parameter in relation to a second other one.
    \n"); cov[1]=1; tj=cptcoveff; @@ -2833,9 +2942,9 @@ void varprob(char optionfilefiname[], do fprintf(ficgp, "**********\n#\n"); - fprintf(fichtm, "\n
    ********** Variable "); + fprintf(fichtmcov, "\n
    ********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - fprintf(fichtm, "**********\n
    "); + fprintf(fichtmcov, "**********\n
    "); fprintf(ficresprobcor, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); @@ -2993,13 +3102,13 @@ void varprob(char optionfilefiname[], do fprintf(ficgp,"\nset parametric;unset label"); fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2); fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65"); - fprintf(fichtm,"\n
    Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\ + fprintf(fichtmcov,"\n
    Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\ :\ %s%d%1d%1d-%1d%1d.png, ",k1,l1,k2,l2,\ subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2,\ subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); - fprintf(fichtm,"\n
    ",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); - fprintf(fichtm,"\n
    Correlation at age %d (%.3f),",(int) age, c12); + fprintf(fichtmcov,"\n
    ",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); + fprintf(fichtmcov,"\n
    Correlation at age %d (%.3f),",(int) age, c12); fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\"",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); @@ -3008,7 +3117,7 @@ void varprob(char optionfilefiname[], do mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); }else{ first=0; - fprintf(fichtm," %d (%.3f),",(int) age, c12); + fprintf(fichtmcov," %d (%.3f),",(int) age, c12); fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\ @@ -3031,7 +3140,8 @@ void varprob(char optionfilefiname[], do fclose(ficresprob); fclose(ficresprobcov); fclose(ficresprobcor); - /* fclose(ficgp);*/ + fflush(ficgp); + fflush(fichtmcov); } @@ -3043,21 +3153,19 @@ void printinghtml(char fileres[], char t double jprev1, double mprev1,double anprev1, \ double jprev2, double mprev2,double anprev2){ int jj1, k1, i1, cpt; - /*char optionfilehtm[FILENAMELENGTH];*/ -/* if((fichtm=fopen(optionfilehtm,"a"))==NULL) { */ -/* printf("Problem with %s \n",optionfilehtm), exit(0); */ -/* fprintf(ficlog,"Problem with %s \n",optionfilehtm), exit(0); */ -/* } */ fprintf(fichtm,"

    • Result files (first order: no variance)

      \n \ - - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): %s
      \n \ - - Estimated transition probabilities over %d (stepm) months: %s
      \n \ - - Stable prevalence in each health state: %s
      \n \ + - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): %s
      \n ", + jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p")); + fprintf(fichtm,"\ + - Estimated transition probabilities over %d (stepm) months: %s
      \n ", + stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij")); + fprintf(fichtm,"\ + - Stable prevalence in each health state: %s
      \n", + subdirf2(fileres,"pl"),subdirf2(fileres,"pl")); + fprintf(fichtm,"\ - Life expectancies by age and initial health status (estepm=%2d months): \ - %s
      \n
      • ", \ - jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p"),\ - stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij"),\ - subdirf2(fileres,"pl"),subdirf2(fileres,"pl"),\ + %s
      \n", estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e")); fprintf(fichtm," \n
      • Graphs

        • "); @@ -3091,28 +3199,32 @@ fprintf(fichtm," \n

        • Graphs fprintf(fichtm,"\n
          - Health life expectancies by age and initial health state (%d): %s%d%d.png
          \ ",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1); } - fprintf(fichtm,"\n
          - Total life expectancy by age and \ -health expectancies in states (1) and (2): %s%d.png
          \ -",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1); } /* end i1 */ }/* End k1 */ fprintf(fichtm,"
        "); - fprintf(fichtm,"\n

      • Result files (second order: variances)

        \n\ - - Parameter file with estimated parameters and covariance matrix: %s
        \n\ - - Variance of one-step probabilities: %s
        \n\ - - Variance-covariance of one-step probabilities: %s
        \n\ - - Correlation matrix of one-step probabilities: %s
        \n\ - - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): %s
        \n\ - - Health expectancies with their variances (no covariance): %s
        \n\ + fprintf(fichtm,"\ +\n

      • Result files (second order: variances)

        \n\ + - Parameter file with estimated parameters and covariance matrix: %s
        \n", rfileres,rfileres); + + fprintf(fichtm," - Variance of one-step probabilities: %s
        \n", + subdirf2(fileres,"prob"),subdirf2(fileres,"prob")); + fprintf(fichtm,"\ + - Variance-covariance of one-step probabilities: %s
        \n", + subdirf2(fileres,"probcov"),subdirf2(fileres,"probcov")); + + fprintf(fichtm,"\ + - Correlation matrix of one-step probabilities: %s
        \n", + subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor")); + fprintf(fichtm,"\ + - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): %s
        \n", + estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v")); + fprintf(fichtm,"\ + - Health expectancies with their variances (no covariance): %s
        \n", + subdirf2(fileres,"t"),subdirf2(fileres,"t")); + fprintf(fichtm,"\ - Standard deviation of stable prevalences: %s
        \n",\ - rfileres,rfileres,\ - subdirf2(fileres,"prob"),subdirf2(fileres,"prob"),\ - subdirf2(fileres,"probcov"),subdirf2(fileres,"probcov"),\ - subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor"),\ - estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"),\ - subdirf2(fileres,"t"),subdirf2(fileres,"t"),\ subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl")); /* if(popforecast==1) fprintf(fichtm,"\n */ @@ -3121,7 +3233,8 @@ health expectancies in states (1) and (2 /*
        ",fileres,fileres,fileres,fileres); */ /* else */ /* fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)

        • \n",popforecast, stepm, model); */ -fprintf(fichtm,"
        • Graphs

          • "); + fflush(fichtm); + fprintf(fichtm,"

          • Graphs

            • "); m=cptcoveff; if (cptcovn < 1) {m=1;ncodemax[1]=1;} @@ -3137,10 +3250,13 @@ fprintf(fichtm,"

            • Graphs"); } for(cpt=1; cpt<=nlstate;cpt++) { - fprintf(fichtm,"
              - Observed and period prevalence (with confident\ -interval) in state (%d): %s%d%d.png
              \ + fprintf(fichtm,"
              - Observed (cross-sectional) and period (incidence based) \ +prevalence (with 95%% confidence interval) in state (%d): %s%d%d.png
              \ ",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1); } + fprintf(fichtm,"\n
              - Total life expectancy by age and \ +health expectancies in states (1) and (2): %s%d.png
              \ +",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1); } /* end i1 */ }/* End k1 */ fprintf(fichtm,"
            "); @@ -3674,9 +3790,14 @@ int fileappend(FILE *fichier, char *opti fflush(fichier); return (1); } + + +/**************** function prwizard **********************/ void prwizard(int ncovmodel, int nlstate, int ndeath, char model[], FILE *ficparo) { + /* Wizard to print covariance matrix template */ + char ca[32], cb[32], cc[32]; int i,j, k, l, li, lj, lk, ll, jj, npar, itimes; int numlinepar; @@ -3790,9 +3911,86 @@ void prwizard(int ncovmodel, int nlstate } /* end k*/ } /*end j */ } /* end i */ - } + } /* 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,"

            • Result files

              \n Force of mortality. Parameters of the Gompertz fit (with confidence interval in brackets):
              "); + fprintf(fichtm," mu(age) =%lf*exp(%lf*(age-%d)) per year

              ",p[1],p[2],agegomp); + for (i=1;i<=2;i++) + fprintf(fichtm," p[%d] = %lf [%f ; %f]
              \n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); + fprintf(fichtm,"

              "); + 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); + +} + + + /***********************************************/ /**************** Main Program *****************/ @@ -3802,8 +4000,12 @@ int main(int argc, char *argv[]) { int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav); int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod; - int jj, imk; + 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 */ /* FILE *ficgp;*/ /*Gnuplot File */ double agedeb, agefin,hf; @@ -3818,7 +4020,7 @@ int main(int argc, char *argv[]) int *indx; char line[MAXLINE], linepar[MAXLINE]; char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; - char pathr[MAXLINE]; + char pathr[MAXLINE], pathimach[MAXLINE]; int firstobs=1, lastobs=10; int sdeb, sfin; /* Status at beginning and end */ int c, h , cpt,l; @@ -3845,9 +4047,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; @@ -3857,15 +4059,16 @@ int main(int argc, char *argv[]) int lstra; long total_usecs; - struct timeval start_time, end_time, curr_time; - struct timezone tzp; - extern int gettimeofday(); - struct tm tmg, tm, *gmtime(), *localtime(); - long time_value; - extern long time(); +/* setlocale (LC_ALL, ""); */ +/* bindtextdomain (PACKAGE, LOCALEDIR); */ +/* textdomain (PACKAGE); */ +/* setlocale (LC_CTYPE, ""); */ +/* setlocale (LC_MESSAGES, ""); */ + /* gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ (void) gettimeofday(&start_time,&tzp); + curr_time=start_time; tm = *localtime(&start_time.tv_sec); tmg = *gmtime(&start_time.tv_sec); strcpy(strstart,asctime(&tm)); @@ -3886,6 +4089,8 @@ int main(int argc, char *argv[]) * printf("tim_value=%d,asctime=%s\n",time_value,strstart); */ + nberr=0; /* Number of errors and warnings */ + nbwarn=0; getcwd(pathcd, size); printf("\n%s\n%s",version,fullversion); @@ -3901,8 +4106,10 @@ int main(int argc, char *argv[]) printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/ /* cutv(path,optionfile,pathtot,'\\');*/ + split(argv[0],pathimach,optionfile,optionfilext,optionfilefiname); + /* strcpy(pathimach,argv[0]); */ split(pathtot,path,optionfile,optionfilext,optionfilefiname); - printf("pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); + printf("pathimach=%s, pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname); chdir(path); strcpy(command,"mkdir "); strcat(command,optionfilefiname); @@ -3927,16 +4134,18 @@ int main(int argc, char *argv[]) } fprintf(ficlog,"Log filename:%s\n",filelog); fprintf(ficlog,"\n%s\n%s",version,fullversion); - fprintf(ficlog,"\nEnter the parameter file name: "); - fprintf(ficlog,"pathtot=%s\n\ + fprintf(ficlog,"\nEnter the parameter file name: \n"); + fprintf(ficlog,"pathimach=%s\npathtot=%s\n\ path=%s \n\ optionfile=%s\n\ optionfilext=%s\n\ - optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); + optionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname); - printf("Localtime (at start):%s",strstart); - fprintf(ficlog,"Localtime (at start): %s",strstart); + printf("Local time (at start):%s",strstart); + fprintf(ficlog,"Local time (at start): %s",strstart); fflush(ficlog); +/* (void) gettimeofday(&curr_time,&tzp); */ +/* printf("Elapsed time %d\n", asc_diff_time(curr_time.tv_sec-start_time.tv_sec,tmpout)); */ /* */ strcpy(fileres,"r"); @@ -3998,51 +4207,144 @@ 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); + else{ + /* 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); + } + 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); - 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); + + 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); + + 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"); } - fprintf(ficparo,"%1d%1d",i1,j1); + } + fflush(ficlog); + + 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 */ + + /* 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); + + matcov=matrix(1,npar,1,npar); + for(i=1; i <=npar; i++){ + fscanf(ficpar,"%s",&str); 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]); + 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(" %lf",param[i][j][k]); - fprintf(ficlog," %lf",param[i][j][k]); + printf(" %.5le",matcov[i][j]); } - else - fprintf(ficlog," %lf",param[i][j][k]); - fprintf(ficparo," %lf",param[i][j][k]); + fprintf(ficlog," %.5le",matcov[i][j]); + fprintf(ficparo," %.5le",matcov[i][j]); } fscanf(ficpar,"\n"); numlinepar++; @@ -4051,115 +4353,33 @@ int main(int argc, char *argv[]) fprintf(ficlog,"\n"); fprintf(ficparo,"\n"); } - } - fflush(ficlog); + 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 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 */ - npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ - - 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); - - 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]); - } - fscanf(ficpar,"\n"); - numlinepar++; - printf("\n"); - fprintf(ficparo,"\n"); - fprintf(ficlog,"\n"); - } - } - fflush(ficlog); - - 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 */ - - /* 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); - - 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]); - } - fscanf(ficpar,"\n"); - numlinepar++; - 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); - - /*-------- data file ----------*/ - if((fic=fopen(datafile,"r"))==NULL) { - printf("Problem with datafile: %s\n", datafile);goto end; - fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end; - } + /*-------- data file ----------*/ + if((fic=fopen(datafile,"r"))==NULL) { + printf("Problem with datafile: %s\n", datafile);goto end; + fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end; + } n= lastobs; severity = vector(1,maxwav); @@ -4182,7 +4402,10 @@ int main(int argc, char *argv[]) i=1; while (fgets(line, MAXLINE, fic) != NULL) { if ((i >= firstobs) && (i <=lastobs)) { - + for(j=0; line[j] != '\n';j++){ /* Untabifies line */ + if(line[j] == '\t') + line[j] = ' '; + } for (j=maxwav;j>=1;j--){ cutv(stra, strb,line,' '); s[j][i]=atoi(strb); strcpy(line,stra); @@ -4328,11 +4551,13 @@ int main(int argc, char *argv[]) s[m][i]=-1; } if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){ + nberr++; printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i); fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i); s[m][i]=-1; } if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){ + nberr++; printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */ @@ -4351,6 +4576,7 @@ int main(int argc, char *argv[]) /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ else { if ((int)andc[i]!=9999){ + nbwarn++; printf("Warning negative age at death: %ld line:%d\n",num[i],i); fprintf(ficlog,"Warning negative age at death: %ld line:%d\n",num[i],i); agev[m][i]=-1; @@ -4387,6 +4613,7 @@ int main(int argc, char *argv[]) for (i=1; i<=imx; i++) { for(m=firstpass; (m<=lastpass); m++){ if (s[m][i] > (nlstate+ndeath)) { + nberr++; printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); goto end; @@ -4401,9 +4628,11 @@ 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); @@ -4457,7 +4686,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); } @@ -4469,12 +4701,26 @@ int main(int argc, char *argv[]) /* fclose(ficgp);*/ /*--------- index.htm --------*/ - strcpy(optionfilehtm,optionfilefiname); + 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); } + strcpy(optionfilehtmcov,optionfilefiname); /* Only for matrix of covariance */ + strcat(optionfilehtmcov,"-cov.htm"); + if((fichtmcov=fopen(optionfilehtmcov,"w"))==NULL) { + printf("Problem with %s \n",optionfilehtmcov), exit(0); + } + else{ + fprintf(fichtmcov,"\nIMaCh Cov %s\n %s
            %s             \ +
            \n\ +Title=%s
            Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
            \n",\ + fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); + } + fprintf(fichtm,"\nIMaCh %s\n %s
            %s             \
            \n\ Title=%s
            Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
            \n\ @@ -4485,18 +4731,14 @@ Title=%s
            Datafile=%s Firstpass=%d La - Log file of the run: %s
            \n\ - Gnuplot file name: %s
            \n\ - Date and time at start: %s
          \n",\ - fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt,\ - model,fileres,fileres,\ + fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\ + fileres,fileres,\ filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart); - /*fclose(fichtm);*/ fflush(fichtm); strcpy(pathr,path); strcat(pathr,optionfilefiname); chdir(optionfilefiname); /* Move to directory named optionfile */ - strcpy(lfileres,fileres); - strcat(lfileres,"/"); - strcat(lfileres,optionfilefiname); /* Calculates basic frequencies. Computes observed prevalence at single age and prints on file fileres'p'. */ @@ -4519,33 +4761,114 @@ Interval (in months) between two waves: 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); - jk=1; - 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) - { + 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 "); + } + + 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); @@ -4559,525 +4882,585 @@ Interval (in months) between two waves: fprintf(ficlog,"\n"); fprintf(ficres,"\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++; + 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("\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"); - for(i=1,k=1;i<=npar;i++){ - /* if (k>nlstate) k=1; - i1=(i-1)/(ncovmodel*nlstate)+1; - fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]); - printf("%s%d%d",alph[k],i1,tab[i]); + + 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 */ - fprintf(ficres,"%3d",i); - if(mle==1) - printf("%3d",i); - fprintf(ficlog,"%3d",i); - for(j=1; j<=i;j++){ - fprintf(ficres," %.5e",matcov[i][j]); - if(mle==1) - printf(" %.5e",matcov[i][j]); - fprintf(ficlog," %.5e",matcov[i][j]); + 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(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); } - fprintf(ficres,"\n"); - if(mle==1) - printf("\n"); - fprintf(ficlog,"\n"); - k++; - } - - 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){ 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) --------------*/ - - 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"); + 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"); + + 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;*/ + 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 */ + agelim=AGESUP; + hstepm=stepsize*YEARM; /* Every year of age */ + hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ - /* hstepm=1; aff par mois*/ + /* 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"); + 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); - - /*---------- Forecasting ------------------*/ - /*if((stepm == 1) && (strcmp(model,".")==0)){*/ - if(prevfcast==1){ - /* if(stepm ==1){*/ + 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){*/ 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(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(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_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------*/ - - 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"); + 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); + + /*------- 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"); - 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); - - /*---------- 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); - - 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); + fclose(ficresvpl); + + /*---------- 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); + + 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_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); - if(erreur >0){ - printf("End of Imach with error or warning %d\n",erreur); - fprintf(ficlog,"End of Imach with error or warning %d\n",erreur); + if((nberr >0) || (nbwarn>0)){ + printf("End of Imach with %d errors and/or %d warnings\n",nberr,nbwarn); + fprintf(ficlog,"End of Imach with %d errors and/or warnings %d\n",nberr,nbwarn); }else{ - printf("End of Imach\n"); - fprintf(ficlog,"End of Imach\n"); + printf("End of Imach\n"); + fprintf(ficlog,"End of Imach\n"); } printf("See log file on %s\n",filelog); - fclose(ficlog); /* gettimeofday(&end_time, (struct timezone*)0);*/ /* after time */ (void) gettimeofday(&end_time,&tzp); tm = *localtime(&end_time.tv_sec); tmg = *gmtime(&end_time.tv_sec); strcpy(strtend,asctime(&tm)); - printf("Localtime at start %s\nLocaltime at end %s",strstart, strtend); - fprintf(ficlog,"Localtime at start %s\nLocal time at end %s",strstart, strtend); - /* printf("Total time used %d Sec\n", asc_time(end_time.tv_sec -start_time.tv_sec);*/ - - printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec); - fprintf(ficlog,"Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec); + printf("Local time at start %s\nLocaltime at end %s",strstart, strtend); + fprintf(ficlog,"Local time at start %s\nLocal time at end %s\n",strstart, strtend); + printf("Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout)); + + printf("Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec); + fprintf(ficlog,"Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout)); + fprintf(ficlog,"Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec); /* printf("Total time was %d uSec.\n", total_usecs);*/ /* if(fileappend(fichtm,optionfilehtm)){ */ fprintf(fichtm,"
          Local time at start %s
          Local time at end %s
          ",strstart, strtend); fclose(fichtm); + fclose(fichtmcov); fclose(ficgp); + fclose(ficlog); /*------ End -----------*/ - end: -#ifdef windows - /* chdir(pathcd);*/ -#endif chdir(path); - /*system("wgnuplot graph.plt");*/ - /*system("../gp37mgw/wgnuplot graph.plt");*/ - /*system("cd ../gp37mgw");*/ - /* system("..\\gp37mgw\\wgnuplot graph.plt");*/ - strcpy(plotcmd,GNUPLOTPROGRAM); + strcpy(plotcmd,"\""); + strcat(plotcmd,pathimach); + strcat(plotcmd,GNUPLOTPROGRAM); + strcat(plotcmd,"\""); strcat(plotcmd," "); strcat(plotcmd,optionfilegnuplot); printf("Starting graphs with: %s",plotcmd);fflush(stdout); - system(plotcmd); + if((outcmd=system(plotcmd)) != 0){ + printf(" Problem with gnuplot\n"); + } printf(" Wait..."); - - /*#ifdef windows*/ while (z[0] != 'q') { /* chdir(path); */ - printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: "); + printf("\nType e to edit output files, g to graph again and q for exiting: "); scanf("%s",z); - if (z[0] == 'c') system("./imach"); - else if (z[0] == 'e') system(optionfilehtm); +/* if (z[0] == 'c') system("./imach"); */ + if (z[0] == 'e') { + printf("Starting browser with: %s",optionfilehtm);fflush(stdout); + system(optionfilehtm); + } else if (z[0] == 'g') system(plotcmd); else if (z[0] == 'q') exit(0); } - /*#endif */ + end: + while (z[0] != 'q') { + printf("\nType q for exiting: "); + scanf("%s",z); + } }