--- imach/src/imach.c 2003/06/25 16:30:45 1.92 +++ imach/src/imach.c 2006/03/15 17:42:26 1.119 @@ -1,6 +1,134 @@ -/* $Id: imach.c,v 1.92 2003/06/25 16:30:45 brouard Exp $ +/* $Id: imach.c,v 1.119 2006/03/15 17:42:26 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.119 2006/03/15 17:42:26 brouard + (Module): Bug if status = -2, the loglikelihood was + computed as likelihood omitting the logarithm. Version O.98e + + Revision 1.118 2006/03/14 18:20:07 brouard + (Module): varevsij Comments added explaining the second + table of variances if popbased=1 . + (Module): Covariances of eij, ekl added, graphs fixed, new html link. + (Module): Function pstamp added + (Module): Version 0.98d + + Revision 1.117 2006/03/14 17:16:22 brouard + (Module): varevsij Comments added explaining the second + table of variances if popbased=1 . + (Module): Covariances of eij, ekl added, graphs fixed, new html link. + (Module): Function pstamp added + (Module): Version 0.98d + + Revision 1.116 2006/03/06 10:29:27 brouard + (Module): Variance-covariance wrong links and + varian-covariance of ej. is needed (Saito). + + Revision 1.115 2006/02/27 12:17:45 brouard + (Module): One freematrix added in mlikeli! 0.98c + + Revision 1.114 2006/02/26 12:57:58 brouard + (Module): Some improvements in processing parameter + filename with strsep. + + Revision 1.113 2006/02/24 14:20:24 brouard + (Module): Memory leaks checks with valgrind and: + datafile was not closed, some imatrix were not freed and on matrix + allocation too. + + Revision 1.112 2006/01/30 09:55:26 brouard + (Module): Back to gnuplot.exe instead of wgnuplot.exe + + Revision 1.111 2006/01/25 20:38:18 brouard + (Module): Lots of cleaning and bugs added (Gompertz) + (Module): Comments can be added in data file. Missing date values + can be a simple dot '.'. + + Revision 1.110 2006/01/25 00:51:50 brouard + (Module): Lots of cleaning and bugs added (Gompertz) + + Revision 1.109 2006/01/24 19:37:15 brouard + (Module): Comments (lines starting with a #) are allowed in data. + + Revision 1.108 2006/01/19 18:05:42 lievre + Gnuplot problem appeared... + To be fixed + + Revision 1.107 2006/01/19 16:20:37 brouard + Test existence of gnuplot in imach path + + Revision 1.106 2006/01/19 13:24:36 brouard + Some cleaning and links added in html output + + Revision 1.105 2006/01/05 20:23:19 lievre + *** empty log message *** + + Revision 1.104 2005/09/30 16:11:43 lievre + (Module): sump fixed, loop imx fixed, and simplifications. + (Module): If the status is missing at the last wave but we know + that the person is alive, then we can code his/her status as -2 + (instead of missing=-1 in earlier versions) and his/her + contributions to the likelihood is 1 - Prob of dying from last + health status (= 1-p13= p11+p12 in the easiest case of somebody in + the healthy state at last known wave). Version is 0.98 + + Revision 1.103 2005/09/30 15:54:49 lievre + (Module): sump fixed, loop imx fixed, and simplifications. + + 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. @@ -100,7 +228,7 @@ hPijx. Also this programme outputs the covariance matrix of the parameters but also - of the life expectancies. It also computes the stable prevalence. + of the life expectancies. It also computes the period (stable) prevalence. Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). Institut national d'études démographiques, Paris. @@ -129,7 +257,7 @@ begin-prev-date,... open gnuplot file open html file - stable prevalence + period (stable) prevalence for age prevalim() h Pij x variance of p varprob @@ -141,7 +269,7 @@ varevsij() if popbased==1 varevsij(,popbased) total life expectancies - Variance of stable prevalence + Variance of period (stable) prevalence end */ @@ -151,18 +279,28 @@ #include #include #include +#include #include -#include +#include +#include +#include +#include +extern int errno; + +/* #include */ #include #include "timeval.h" +/* #include */ +/* #define _(String) gettext (String) */ + #define MAXLINE 256 + #define GNUPLOTPROGRAM "gnuplot" /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ #define FILENAMELENGTH 132 -/*#define DEBUG*/ -/*#define windows*/ + #define GLOCK_ERROR_NOPATH -1 /* empty path */ #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ @@ -177,19 +315,24 @@ #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 CHARSEPARATOR "/" #define ODIRSEPARATOR '\\' #else #define DIRSEPARATOR '\\' +#define CHARSEPARATOR "\\" #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.92 2003/06/25 16:30:45 brouard Exp $ */ +/* $Id: imach.c,v 1.119 2006/03/15 17:42:26 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.96a, June 2003, INED-EUROREVES "; -char fullversion[]="$Revision: 1.92 $ $Date: 2003/06/25 16:30:45 $"; +char version[]="Imach version 0.98e, March 2006, INED-EUROREVES-Institut de longevite "; +char fullversion[]="$Revision: 1.119 $ $Date: 2006/03/15 17:42:26 $"; +char strstart[80]; +char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ int nvar; int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; @@ -202,6 +345,7 @@ int popbased=0; int *wav; /* Number of waves for this individuual 0 is possible */ int maxwav; /* Maxim number of waves */ int jmin, jmax; /* min, max spacing between 2 waves */ +int ijmin, ijmax; /* Individuals having jmin and jmax */ int gipmx, gsw; /* Global variables on the number of contributions to the likelihood and the sum of weights (done by funcone)*/ int mle, weightopt; @@ -212,12 +356,13 @@ int **bh; /* bh[mi][i] is the bias (+ or double jmean; /* Mean space between 2 waves */ double **oldm, **newm, **savm; /* Working pointers to matrices */ double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ -FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; +FILE *fic,*ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; FILE *ficlog, *ficrespow; int globpr; /* Global variable for printing or not */ double fretone; /* Only one call to likelihood */ long ipmx; /* Number of contributions */ double sw; /* Sum of weights */ +char filerespow[FILENAMELENGTH]; char fileresilk[FILENAMELENGTH]; /* File of individual contributions to the likelihood */ FILE *ficresilk; FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor; @@ -225,19 +370,23 @@ FILE *ficresprobmorprev; FILE *fichtm, *fichtmcov; /* Html File */ FILE *ficreseij; char filerese[FILENAMELENGTH]; +FILE *ficresstdeij; +char fileresstde[FILENAMELENGTH]; +FILE *ficrescveij; +char filerescve[FILENAMELENGTH]; FILE *ficresvij; char fileresv[FILENAMELENGTH]; FILE *ficresvpl; char fileresvpl[FILENAMELENGTH]; char title[MAXLINE]; char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH]; -char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH]; -char tmpout[FILENAMELENGTH]; +char plotcmd[FILENAMELENGTH], pplotcmd[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]; @@ -253,6 +402,9 @@ long time_value; extern long time(); char strcurr[80], strfor[80]; +char *endptr; +long lval; + #define NR_END 1 #define FREE_ARG char* #define FTOL 1.0e-10 @@ -280,9 +432,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; @@ -290,15 +443,17 @@ 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; int **s; /* Status */ double *agedc, **covar, idx; int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff; +double *lsurv, *lpop, *tpop; double ftol=FTOL; /* Tolerance for computing Max Likelihood */ double ftolhess; /* Tolerance for computing hessian */ @@ -306,13 +461,17 @@ 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 */ l1 = strlen(path ); /* length of path */ if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH ); ss= strrchr( path, DIRSEPARATOR ); /* find last / */ - if ( ss == NULL ) { /* no directory, so use current */ + if ( ss == NULL ) { /* no directory, so determine current directory */ + strcpy( name, path ); /* we got the fullname name because no directory */ /*if(strrchr(path, ODIRSEPARATOR )==NULL) printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/ /* get current working directory */ @@ -320,7 +479,8 @@ static int split( char *path, char *dirc if ( getcwd( dirc, FILENAME_MAX ) == NULL ) { return( GLOCK_ERROR_GETCWD ); } - strcpy( name, path ); /* we've got it */ + /* got dirc from getcwd*/ + printf(" DIRC = %s \n",dirc); } else { /* strip direcotry from path */ ss++; /* after this, the filename */ l2 = strlen( ss ); /* length of filename */ @@ -328,21 +488,25 @@ static int split( char *path, char *dirc strcpy( name, ss ); /* save file name */ strncpy( dirc, path, l1 - l2 ); /* now the directory */ dirc[l1-l2] = 0; /* add zero */ + printf(" DIRC2 = %s \n",dirc); } + /* We add a separator at the end of dirc if not exists */ l1 = strlen( dirc ); /* length of directory */ - /*#ifdef windows - if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; } -#else - if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; } -#endif - */ + if( dirc[l1-1] != DIRSEPARATOR ){ + dirc[l1] = DIRSEPARATOR; + dirc[l1+1] = 0; + printf(" DIRC3 = %s \n",dirc); + } 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 */ } @@ -375,8 +539,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; @@ -557,6 +721,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; @@ -770,9 +969,10 @@ void powell(double p[], double **xi, int last_time=curr_time; (void) gettimeofday(&curr_time,&tzp); printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout); - fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); + /* fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); - for (i=1;i<=n;i++) { + */ + for (i=1;i<=n;i++) { printf(" %d %.12f",i, p[i]); fprintf(ficlog," %d %.12lf",i, p[i]); fprintf(ficrespow," %.12lf", p[i]); @@ -782,11 +982,11 @@ void powell(double p[], double **xi, int fprintf(ficrespow,"\n");fflush(ficrespow); if(*iter <=3){ tm = *localtime(&curr_time.tv_sec); - strcpy(strcurr,asctime(&tmf)); + strcpy(strcurr,asctime(&tm)); /* asctime_r(&tm,strcurr); */ - forecast_time=curr_time; + forecast_time=curr_time; itmp = strlen(strcurr); - if(strcurr[itmp-1]=='\n') + 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); @@ -798,8 +998,8 @@ void powell(double p[], double **xi, int 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 or\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 or\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); + 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++) { @@ -893,7 +1093,7 @@ void powell(double p[], double **xi, int } } -/**** Prevalence limit (stable prevalence) ****************/ +/**** Prevalence limit (stable or period prevalence) ****************/ double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij) { @@ -964,57 +1164,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. @@ -1153,15 +1353,16 @@ 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 . + /* 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 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 @@ -1181,7 +1382,28 @@ double func( double *x) lower mortality. */ lli=log(out[s1][s2] - savm[s1][s2]); - }else{ + + + } else if (s2==-2) { + for (j=1,survp=0. ; j<=nlstate; j++) + survp += out[s1][j]; + lli= log(survp); + } + +/* else if (s2==-4) { */ +/* for (j=3,survp=0. ; j<=nlstate; j++) */ +/* survp += out[s1][j]; */ +/* lli= survp; */ +/* } */ + +/* else if (s2==-5) { */ +/* for (j=1,survp=0. ; j<=2; j++) */ +/* survp += out[s1][j]; */ +/* lli= survp; */ +/* } */ + + + else{ lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ /* 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 */ } @@ -1214,31 +1436,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; @@ -1265,30 +1466,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; @@ -1416,7 +1597,11 @@ double funcone( double *x) */ if( s2 > nlstate && (mle <5) ){ /* Jackson */ lli=log(out[s1][s2] - savm[s1][s2]); - } else if (mle==1){ + } else if (s2==-2) { + for (j=1,survp=0. ; j<=nlstate; j++) + survp += out[s1][j]; + lli= log(survp); + }else if (mle==1){ lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ } else if(mle==2){ 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 */ @@ -1432,8 +1617,8 @@ double funcone( double *x) ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; /* printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ if(globpr){ - fprintf(ficresilk,"%9d %6d %1d %1d %1d %1d %3d %10.6f %6.4f\ - %10.6f %10.6f %10.6f ", \ + fprintf(ficresilk,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\ + %11.6f %11.6f %11.6f ", \ num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i], 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); for(k=1,llt=0.,l=0.; k<=nlstate; k++){ @@ -1454,35 +1639,8 @@ double funcone( double *x) return -l; } -char *subdirf(char fileres[]) -{ - /* Caution optionfilefiname is hidden */ - 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 @@ -1525,7 +1683,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++) @@ -1545,6 +1703,7 @@ void mlikeli(FILE *ficres,double p[], in powell(p,xi,npar,ftol,&iter,&fret,func); + free_matrix(xi,1,npar,1,npar); fclose(ficrespow); printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p)); fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); @@ -1560,11 +1719,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"); @@ -1572,9 +1731,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++) { @@ -1582,7 +1743,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]);*/ } @@ -1653,14 +1815,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; @@ -1700,7 +1862,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; @@ -1809,8 +1971,13 @@ void lubksb(double **a, int n, int *indx } } +void pstamp(FILE *fichier) +{ + fprintf(fichier,"# %s.%s\n#%s\n#%s\n# %s", optionfilefiname,optionfilext,version,fullversion,strstart); +} + /************ Frequencies ********************/ -void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint) +void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[]) { /* Some frequencies */ int i, m, jk, k1,i1, j1, bool, z1,z2,j; @@ -1818,7 +1985,6 @@ void freqsummary(char fileres[], int ia double ***freq; /* Frequencies */ double *pp, **prop; double pos,posprop, k2, dateintsum=0,k2cpt=0; - FILE *ficresp; char fileresp[FILENAMELENGTH]; pp=vector(1,nlstate); @@ -1830,7 +1996,7 @@ void freqsummary(char fileres[], int ia fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); exit(0); } - freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3); + freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3); j1=0; j=cptcoveff; @@ -1843,8 +2009,8 @@ void freqsummary(char fileres[], int ia j1++; /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); scanf("%d", i);*/ - for (i=-1; i<=nlstate+ndeath; i++) - for (jk=-1; jk<=nlstate+ndeath; jk++) + for (i=-5; i<=nlstate+ndeath; i++) + for (jk=-5; jk<=nlstate+ndeath; jk++) for(m=iagemin; m <= iagemax+3; m++) freq[i][jk][m]=0; @@ -1883,7 +2049,7 @@ void freqsummary(char fileres[], int ia } /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ - + pstamp(ficresp); if (cptcovn>0) { fprintf(ficresp, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); @@ -1969,7 +2135,7 @@ void freqsummary(char fileres[], int ia dateintmean=dateintsum/k2cpt; fclose(ficresp); - free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3); + free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3); free_vector(pp,1,nlstate); free_matrix(prop,1,nlstate,iagemin, iagemax+3); /* End of Freq */ @@ -2078,7 +2244,7 @@ void concatwav(int wav[], int **dh, int mi=0; m=firstpass; while(s[m][i] <= nlstate){ - if(s[m][i]>=1) + if(s[m][i]>=1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5) mw[++mi][i]=m; if(m >=lastpass) break; @@ -2096,11 +2262,11 @@ void concatwav(int wav[], int **dh, int 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); + printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i); first=1; } if(first==1){ - fprintf(ficlog,"Warning! None valid information for:%ld line=%d (skipped)\n",num[i],i); + fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i); } } /* end mi==0 */ } /* End individuals */ @@ -2118,13 +2284,19 @@ void concatwav(int wav[], int **dh, int 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; /* 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(" 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. You MUST fix the contradiction between dates.\n",stepm); 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); + 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. You MUST fix the contradiction between dates.\n",stepm); } k=k+1; - if (j >= jmax) jmax=j; - if (j <= jmin) jmin=j; + if (j >= jmax){ + jmax=j; + ijmax=i; + } + if (j <= jmin){ + jmin=j; + ijmin=i; + } sum=sum+j; /*if (j<0) printf("j=%d num=%d \n",j,i);*/ /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ @@ -2132,10 +2304,17 @@ void concatwav(int wav[], int **dh, int } else{ j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); - /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ +/* if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */ + k=k+1; - if (j >= jmax) jmax=j; - else if (j <= jmin)jmin=j; + if (j >= jmax) { + jmax=j; + ijmax=i; + } + else if (j <= jmin){ + jmin=j; + ijmin=i; + } /* 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){ @@ -2178,8 +2357,8 @@ void concatwav(int wav[], int **dh, int } /* end wave */ } jmean=sum/k; - printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean); - fprintf(ficlog,"Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean); + printf("Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, num[ijmin], jmax, num[ijmax], jmean); + fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean); } /*********** Tricode ****************************/ @@ -2225,7 +2404,7 @@ void tricode(int *Tvar, int **nbcode, in for (k=0; k< maxncov; k++) Ndum[k]=0; for (i=1; i<=ncovmodel-2; i++) { - /* Listing of all covariables in staement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ + /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ ij=Tvar[i]; Ndum[ij]++; } @@ -2243,31 +2422,151 @@ void tricode(int *Tvar, int **nbcode, in /*********** Health Expectancies ****************/ -void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm,double delti[],double **matcov ) +void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] ) + +{ + /* Health expectancies, no variances */ + int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2; + double age, agelim, hf; + double ***p3mat; + double eip; + + pstamp(ficreseij); + fprintf(ficreseij,"# (a) Life expectancies by health status at initial age and (b) health expectancies by health status at initial age\n"); + fprintf(ficreseij,"# Age"); + for(i=1; i<=nlstate;i++){ + for(j=1; j<=nlstate;j++){ + fprintf(ficreseij," e%1d%1d ",i,j); + } + fprintf(ficreseij," e%1d. ",i); + } + fprintf(ficreseij,"\n"); + + + if(estepm < stepm){ + printf ("Problem %d lower than %d\n",estepm, stepm); + } + else hstepm=estepm; + /* We compute the life expectancy from trapezoids spaced every estepm months + * This is mainly to measure the difference between two models: for example + * if stepm=24 months pijx are given only every 2 years and by summing them + * we are calculating an estimate of the Life Expectancy assuming a linear + * progression in between and thus overestimating or underestimating according + * to the curvature of the survival function. If, for the same date, we + * estimate the model with stepm=1 month, we can keep estepm to 24 months + * to compare the new estimate of Life expectancy with the same linear + * hypothesis. A more precise result, taking into account a more precise + * curvature will be obtained if estepm is as small as stepm. */ + + /* For example we decided to compute the life expectancy with the smallest unit */ + /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. + nhstepm is the number of hstepm from age to agelim + nstepm is the number of stepm from age to agelin. + Look at hpijx to understand the reason of that which relies in memory size + and note for a fixed period like estepm months */ + /* We decided (b) to get a life expectancy respecting the most precise curvature of the + survival function given by stepm (the optimization length). Unfortunately it + means that if the survival funtion is printed only each two years of age and if + you sum them up and add 1 year (area under the trapezoids) you won't get the same + results. So we changed our mind and took the option of the best precision. + */ + hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ + + agelim=AGESUP; + /* nhstepm age range expressed in number of stepm */ + nstepm=(int) rint((agelim-age)*YEARM/stepm); + /* Typically if 20 years nstepm = 20*12/6=40 stepm */ + /* if (stepm >= YEARM) hstepm=1;*/ + nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ + p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + + for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ + /* Computed by stepm unit matrices, product of hstepm matrices, stored + in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ + + hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); + + hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ + + printf("%d|",(int)age);fflush(stdout); + fprintf(ficlog,"%d|",(int)age);fflush(ficlog); + + /* Computing expectancies */ + for(i=1; i<=nlstate;i++) + for(j=1; j<=nlstate;j++) + for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ + eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf; + + /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ + + } + + fprintf(ficreseij,"%3.0f",age ); + for(i=1; i<=nlstate;i++){ + eip=0; + for(j=1; j<=nlstate;j++){ + eip +=eij[i][j][(int)age]; + fprintf(ficreseij,"%9.4f", eij[i][j][(int)age] ); + } + fprintf(ficreseij,"%9.4f", eip ); + } + fprintf(ficreseij,"\n"); + + } + free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + printf("\n"); + fprintf(ficlog,"\n"); + +} + +void cvevsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] ) { - /* Health expectancies */ - int i, j, nhstepm, hstepm, h, nstepm, k, cptj; + /* Covariances of health expectancies eij and of total life expectancies according + to initial status i, ei. . + */ + int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji; double age, agelim, hf; - double ***p3mat,***varhe; + double ***p3matp, ***p3matm, ***varhe; double **dnewm,**doldm; - double *xp; + double *xp, *xm; double **gp, **gm; double ***gradg, ***trgradg; int theta; + double eip, vip; + varhe=ma3x(1,nlstate*nlstate,1,nlstate*nlstate,(int) bage, (int) fage); xp=vector(1,npar); + xm=vector(1,npar); dnewm=matrix(1,nlstate*nlstate,1,npar); doldm=matrix(1,nlstate*nlstate,1,nlstate*nlstate); - fprintf(ficreseij,"# Health expectancies\n"); - fprintf(ficreseij,"# Age"); - for(i=1; i<=nlstate;i++) + pstamp(ficresstdeij); + fprintf(ficresstdeij,"# Health expectancies with standard errors\n"); + fprintf(ficresstdeij,"# Age"); + for(i=1; i<=nlstate;i++){ for(j=1; j<=nlstate;j++) - fprintf(ficreseij," %1d-%1d (SE)",i,j); - fprintf(ficreseij,"\n"); + fprintf(ficresstdeij," e%1d%1d (SE)",i,j); + fprintf(ficresstdeij," e%1d. ",i); + } + fprintf(ficresstdeij,"\n"); + pstamp(ficrescveij); + fprintf(ficrescveij,"# Subdiagonal matrix of covariances of health expectancies by age: cov(eij,ekl)\n"); + fprintf(ficrescveij,"# Age"); + for(i=1; i<=nlstate;i++) + for(j=1; j<=nlstate;j++){ + cptj= (j-1)*nlstate+i; + for(i2=1; i2<=nlstate;i2++) + for(j2=1; j2<=nlstate;j2++){ + cptj2= (j2-1)*nlstate+i2; + if(cptj2 <= cptj) + fprintf(ficrescveij," %1d%1d,%1d%1d",i,j,i2,j2); + } + } + fprintf(ficrescveij,"\n"); + if(estepm < stepm){ printf ("Problem %d lower than %d\n",estepm, stepm); } @@ -2297,77 +2596,64 @@ void evsij(char fileres[], double ***eij */ hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ + /* If stepm=6 months */ + /* nhstepm age range expressed in number of stepm */ agelim=AGESUP; - for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ - /* nhstepm age range expressed in number of stepm */ - nstepm=(int) rint((agelim-age)*YEARM/stepm); - /* Typically if 20 years nstepm = 20*12/6=40 stepm */ - /* if (stepm >= YEARM) hstepm=1;*/ - nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ - p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - gradg=ma3x(0,nhstepm,1,npar,1,nlstate*nlstate); - gp=matrix(0,nhstepm,1,nlstate*nlstate); - gm=matrix(0,nhstepm,1,nlstate*nlstate); + nstepm=(int) rint((agelim-age)*YEARM/stepm); + /* Typically if 20 years nstepm = 20*12/6=40 stepm */ + /* if (stepm >= YEARM) hstepm=1;*/ + nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ + + p3matp=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + p3matm=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + gradg=ma3x(0,nhstepm,1,npar,1,nlstate*nlstate); + trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar); + gp=matrix(0,nhstepm,1,nlstate*nlstate); + gm=matrix(0,nhstepm,1,nlstate*nlstate); + + for (age=bage; age<=fage; age ++){ /* Computed by stepm unit matrices, product of hstepm matrices, stored in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ - hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij); - hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ - /* Computing Variances of health expectancies */ - - for(theta=1; theta <=npar; theta++){ + /* Computing Variances of health expectancies */ + /* Gradient is computed with plus gp and minus gm. Code is duplicated in order to + decrease memory allocation */ + for(theta=1; theta <=npar; theta++){ for(i=1; i<=npar; i++){ xp[i] = x[i] + (i==theta ?delti[theta]:0); + xm[i] = x[i] - (i==theta ?delti[theta]:0); } - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); + hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij); + hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij); - cptj=0; for(j=1; j<= nlstate; j++){ for(i=1; i<=nlstate; i++){ - cptj=cptj+1; - for(h=0, gp[h][cptj]=0.; h<=nhstepm-1; h++){ - gp[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.; + for(h=0; h<=nhstepm-1; h++){ + gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.; + gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.; } } } - - for(i=1; i<=npar; i++) - xp[i] = x[i] - (i==theta ?delti[theta]:0); - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); - - cptj=0; - for(j=1; j<= nlstate; j++){ - for(i=1;i<=nlstate;i++){ - cptj=cptj+1; - for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){ - - gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.; - } - } - } - for(j=1; j<= nlstate*nlstate; j++) + for(ij=1; ij<= nlstate*nlstate; ij++) for(h=0; h<=nhstepm-1; h++){ - gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; + gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta]; } - } - -/* End theta */ - - trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar); - - for(h=0; h<=nhstepm-1; h++) + }/* End theta */ + + + for(h=0; h<=nhstepm-1; h++) for(j=1; j<=nlstate*nlstate;j++) for(theta=1; theta <=npar; theta++) trgradg[h][j][theta]=gradg[h][theta][j]; - + - for(i=1;i<=nlstate*nlstate;i++) - for(j=1;j<=nlstate*nlstate;j++) - varhe[i][j][(int)age] =0.; + for(ij=1;ij<=nlstate*nlstate;ij++) + for(ji=1;ji<=nlstate*nlstate;ji++) + varhe[ij][ji][(int)age] =0.; printf("%d|",(int)age);fflush(stdout); fprintf(ficlog,"%d|",(int)age);fflush(ficlog); @@ -2375,39 +2661,60 @@ void evsij(char fileres[], double ***eij for(k=0;k<=nhstepm-1;k++){ matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov); matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]); - for(i=1;i<=nlstate*nlstate;i++) - for(j=1;j<=nlstate*nlstate;j++) - varhe[i][j][(int)age] += doldm[i][j]*hf*hf; + for(ij=1;ij<=nlstate*nlstate;ij++) + for(ji=1;ji<=nlstate*nlstate;ji++) + varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf; } } /* Computing expectancies */ + hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ - eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf; + eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[i][j][h+1])/2.0*hf; -/* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ + /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ } - fprintf(ficreseij,"%3.0f",age ); - cptj=0; + fprintf(ficresstdeij,"%3.0f",age ); + for(i=1; i<=nlstate;i++){ + eip=0.; + vip=0.; + for(j=1; j<=nlstate;j++){ + eip += eij[i][j][(int)age]; + for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */ + vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age]; + fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) ); + } + fprintf(ficresstdeij," %9.4f (%.4f)", eip, sqrt(vip)); + } + fprintf(ficresstdeij,"\n"); + + fprintf(ficrescveij,"%3.0f",age ); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++){ - cptj++; - fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) ); + cptj= (j-1)*nlstate+i; + for(i2=1; i2<=nlstate;i2++) + for(j2=1; j2<=nlstate;j2++){ + cptj2= (j2-1)*nlstate+i2; + if(cptj2 <= cptj) + fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]); + } } - fprintf(ficreseij,"\n"); + fprintf(ficrescveij,"\n"); - free_matrix(gm,0,nhstepm,1,nlstate*nlstate); - free_matrix(gp,0,nhstepm,1,nlstate*nlstate); - free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*nlstate); - free_ma3x(trgradg,0,nhstepm,1,nlstate*nlstate,1,npar); - free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); } + free_matrix(gm,0,nhstepm,1,nlstate*nlstate); + free_matrix(gp,0,nhstepm,1,nlstate*nlstate); + free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*nlstate); + free_ma3x(trgradg,0,nhstepm,1,nlstate*nlstate,1,npar); + free_ma3x(p3matm,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + free_ma3x(p3matp,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); printf("\n"); fprintf(ficlog,"\n"); + free_vector(xm,1,npar); free_vector(xp,1,npar); free_matrix(dnewm,1,nlstate*nlstate,1,npar); free_matrix(doldm,1,nlstate*nlstate,1,nlstate*nlstate); @@ -2415,7 +2722,7 @@ void evsij(char fileres[], double ***eij } /************ Variance ******************/ -void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav) +void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[]) { /* Variance of health expectancies */ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ @@ -2466,7 +2773,9 @@ void varevsij(char optionfilefiname[], d fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev); } printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); + fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); + pstamp(ficresprobmorprev); fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm); fprintf(ficresprobmorprev,"# Age cov=%-d",ij); for(j=nlstate+1; j<=(nlstate+ndeath);j++){ @@ -2476,16 +2785,21 @@ void varevsij(char optionfilefiname[], d } fprintf(ficresprobmorprev,"\n"); fprintf(ficgp,"\n# Routine varevsij"); + /* fprintf(fichtm, "#Local time at start: %s", strstart);*/ fprintf(fichtm,"\n

  • Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)

    • \n"); fprintf(fichtm,"\n
    %s
    \n",digitp); /* } */ varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); - - fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are the stable prevalence in health states i\n"); + pstamp(ficresvij); + fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are "); + if(popbased==1) + fprintf(ficresvij,"the age specific prevalence observed in the population i.e cross-sectionally\n in each health state (popbased=1)"); + else + fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n"); fprintf(ficresvij,"# Age"); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) - fprintf(ficresvij," Cov(e%1d, e%1d)",i,j); + fprintf(ficresvij," Cov(e.%1d, e.%1d)",i,j); fprintf(ficresvij,"\n"); xp=vector(1,npar); @@ -2713,7 +3027,7 @@ void varevsij(char optionfilefiname[], d } /* end varevsij */ /************ Variance of prevlim ******************/ -void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij) +void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, char strstart[]) { /* Variance of prevalence limit */ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ @@ -2726,8 +3040,9 @@ void varprevlim(char fileres[], double * double **gradg, **trgradg; double age,agelim; int theta; - - fprintf(ficresvpl,"# Standard deviation of stable prevalences \n"); + + pstamp(ficresvpl); + fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n"); fprintf(ficresvpl,"# Age"); for(i=1; i<=nlstate;i++) fprintf(ficresvpl," %1d-%1d",i,i); @@ -2796,7 +3111,7 @@ void varprevlim(char fileres[], double * } /************ Variance of one-step probabilities ******************/ -void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax) +void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax, char strstart[]) { int i, j=0, i1, k1, l1, t, tj; int k2, l2, j1, z1; @@ -2841,13 +3156,15 @@ void varprob(char optionfilefiname[], do fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); - + pstamp(ficresprob); fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n"); fprintf(ficresprob,"# Age"); + pstamp(ficresprobcov); fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); fprintf(ficresprobcov,"# Age"); + pstamp(ficresprobcor); fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); - fprintf(ficresprobcov,"# Age"); + fprintf(ficresprobcor,"# Age"); for(i=1; i<=nlstate;i++) @@ -2870,8 +3187,8 @@ 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",optionfilehtmcov); - fprintf(fichtmcov,"\n

    Computing matrix of variance-covariance of step probabilities

    \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.\ @@ -2904,9 +3221,9 @@ To be simple, these graphs help to under 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]]); @@ -3098,6 +3415,8 @@ To be simple, these graphs help to under } free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); + free_matrix(doldm,1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); + free_matrix(dnewm,1,(nlstate)*(nlstate+ndeath),1,npar); free_vector(xp,1,npar); fclose(ficresprob); fclose(ficresprobcov); @@ -3115,23 +3434,25 @@ 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 \ - - 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"),\ + fprintf(fichtm,""); + 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 ", + 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,"\ + - Period (stable) prevalence in each health state: %s
        \n", + subdirf2(fileres,"pl"),subdirf2(fileres,"pl")); + fprintf(fichtm,"\ + - (a) Life expectancies by health status at initial age, (b) health expectancies by health status at initial age: ei., eij (estepm=%2d months): \ + %s
        \n
        • ", estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e")); + fprintf(fichtm," \n
        • Graphs

          • "); m=cptcoveff; @@ -3148,43 +3469,55 @@ fprintf(fichtm," \n

          • Graphs fprintf(fichtm," ************\n
            "); } /* Pij */ - fprintf(fichtm,"
            - Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s%d1.png
            \ -",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); + fprintf(fichtm,"
            - Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s%d1.png
            \ +",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); /* Quasi-incidences */ fprintf(fichtm,"
            - Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\ - before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: %s%d2.png
            \ -",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); - /* Stable prevalence in each health state */ + before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: %s%d2.png
            \ +",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); + /* Period (stable) prevalence in each health state */ for(cpt=1; cpt- Stable prevalence in each health state : p%s%d%d.png
            \ -",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1); + fprintf(fichtm,"
            - Period (stable) prevalence in each health state : %s%d%d.png
            \ +",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1); } for(cpt=1; cpt<=nlstate;cpt++) { - 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
            - Life expectancy by health state (%d) at initial age and its decomposition into health expectancies : %s%d%d.png
            \ +",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,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\ - - 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"),\ + 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 health expectancies by age and initial health status (cov(eij,ekl)(estepm=%2d months): \ + %s
          \n
          • ", + estepm,subdirf2(fileres,"cve"),subdirf2(fileres,"cve")); + fprintf(fichtm,"\ + - (a) Health expectancies by health status at initial age (eij) and standard errors (in parentheses) (b) life expectancies and standard errors (ei.=ei1+ei2+...)(estepm=%2d months): \ + %s
          \n", + estepm,subdirf2(fileres,"stde"),subdirf2(fileres,"stde")); + fprintf(fichtm,"\ + - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences (i.e population based) (estepm=%d months): %s
          \n", + estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v")); + fprintf(fichtm,"\ + - Total life expectancy and total health expectancies to be spent in each health state e.j with their standard errors: %s
          \n", + subdirf2(fileres,"t"),subdirf2(fileres,"t")); + fprintf(fichtm,"\ + - Standard deviation of period (stable) prevalences: %s
          \n",\ subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl")); /* if(popforecast==1) fprintf(fichtm,"\n */ @@ -3193,7 +3526,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;} @@ -3209,10 +3543,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,"
              "); @@ -3252,7 +3589,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); else fprintf(ficgp," \%%*lf (\%%*lf)"); } - fprintf(ficgp,"\" t\"Stable prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); + fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); else fprintf(ficgp," \%%*lf (\%%*lf)"); @@ -3300,7 +3637,8 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u for (k1=1; k1<= m ; k1 ++) { for (cpt=1; cpt<= nlstate ; cpt ++) { - k=2+nlstate*(2*cpt-2); + /* k=2+nlstate*(2*cpt-2); */ + k=2+(nlstate+1)*(cpt-1); fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1); fprintf(ficgp,"set ter png small\n\ set size 0.65,0.65\n\ @@ -3314,9 +3652,11 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u */ for (i=1; i< nlstate ; i ++) { - fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1); + fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+i,cpt,i+1); + /* fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/ } + fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+nlstate,cpt); } } @@ -3746,9 +4086,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; @@ -3862,9 +4207,95 @@ 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=0; 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=1;i<=imx ; i++) + { + if (cens[i] == 1 && wav[i]>1) + A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); + + if (cens[i] == 0 && wav[i]>1) + A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))) + +log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM); + + /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */ + if (wav[i] > 1 ) { /* ??? */ + 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,double agemortsup){ + int i,k; + + 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,"
              "); + +fprintf(fichtm,"

              • Life table

                \n
                "); + + fprintf(fichtm,"\nAge lx qx d(x,x+1) Lx Tx e
                "); + + for (k=agegomp;k<(agemortsup-2);k++) + fprintf(fichtm,"%d %.0lf %lf %.0lf %.0lf %.0lf %lf
                \n",k,lsurv[k],p[1]*exp(p[2]*(k-agegomp)),(p[1]*exp(p[2]*(k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]); + + + 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); + +} + + + /***********************************************/ @@ -3875,10 +4306,17 @@ 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 linei, month, year,iout; + 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]; + char dummy[]=" "; /* FILE *fichtm; *//* Html File */ /* FILE *ficgp;*/ /*Gnuplot File */ + struct stat info; double agedeb, agefin,hf; double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; @@ -3891,7 +4329,8 @@ 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]; + char **bp, *tok, *val; /* pathtot */ int firstobs=1, lastobs=10; int sdeb, sfin; /* Status at beginning and end */ int c, h , cpt,l; @@ -3901,6 +4340,8 @@ int main(int argc, char *argv[]) int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ int mobilav=0,popforecast=0; int hstepm, nhstepm; + int agemortsup; + float sumlpop=0.; double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000; double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000; @@ -3918,19 +4359,25 @@ 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; char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80]; - char strstart[80], *strt, strtend[80]; + char *strt, strtend[80]; char *stratrunc; int lstra; long total_usecs; +/* 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; @@ -3961,7 +4408,17 @@ int main(int argc, char *argv[]) printf("\n%s\n%s",version,fullversion); if(argc <=1){ printf("\nEnter the parameter file name: "); - scanf("%s",pathtot); + fgets(pathr,FILENAMELENGTH,stdin); + i=strlen(pathr); + if(pathr[i-1]=='\n') + pathr[i-1]='\0'; + for (tok = pathr; tok != NULL; ){ + printf("Pathr |%s|\n",pathr); + while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0'); + printf("val= |%s| pathr=%s\n",val,pathr); + strcpy (pathtot, val); + if(pathr[0] == '\0') break; /* Dirty */ + } } else{ strcpy(pathtot,argv[1]); @@ -3971,9 +4428,17 @@ int main(int argc, char *argv[]) printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/ /* cutv(path,optionfile,pathtot,'\\');*/ + /* Split argv[0], imach program to get pathimach */ + printf("\nargv[0]=%s argv[1]=%s, \n",argv[0],argv[1]); + split(argv[0],pathimach,optionfile,optionfilext,optionfilefiname); + printf("\nargv[0]=%s pathimach=%s, \noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",argv[0],pathimach,optionfile,optionfilext,optionfilefiname); + /* strcpy(pathimach,argv[0]); */ + /* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */ split(pathtot,path,optionfile,optionfilext,optionfilefiname); - printf("pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); - chdir(path); + printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); + chdir(path); /* Can be a relative path */ + if(getcwd(pathcd,MAXLINE) > 0) /* So pathcd is the full path */ + printf("Current directory %s!\n",pathcd); strcpy(command,"mkdir "); strcat(command,optionfilefiname); if((outcmd=system(command)) != 0){ @@ -3997,15 +4462,15 @@ 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)); */ @@ -4070,167 +4535,177 @@ int main(int argc, char *argv[]) ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ - + npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ + + delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); + delti=delti3[1][1]; + /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ if(mle==-1){ /* Print a wizard for help writing covariance matrix */ prwizard(ncovmodel, nlstate, ndeath, model, ficparo); printf(" You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso); fprintf(ficlog," You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso); + free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); fclose (ficparo); fclose (ficlog); exit(0); } - /* Read guess parameters */ - /* Reads comments: lines beginning with '#' */ - while((c=getc(ficpar))=='#' && c!= EOF){ - ungetc(c,ficpar); - fgets(line, MAXLINE, ficpar); - numlinepar++; - puts(line); - fputs(line,ficparo); - fputs(line,ficlog); + else if(mle==-3) { + prwizard(ncovmodel, nlstate, ndeath, model, ficparo); + printf(" You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso); + fprintf(ficlog," You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso); + param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); + matcov=matrix(1,npar,1,npar); } - ungetc(c,ficpar); - - param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); - for(i=1; i <=nlstate; i++){ - j=0; - for(jj=1; jj <=nlstate+ndeath; jj++){ - if(jj==i) continue; - j++; - fscanf(ficpar,"%1d%1d",&i1,&j1); - if ((i1 != i) && (j1 != j)){ - printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); - exit(1); - } - fprintf(ficparo,"%1d%1d",i1,j1); - if(mle==1) - printf("%1d%1d",i,j); - fprintf(ficlog,"%1d%1d",i,j); - for(k=1; k<=ncovmodel;k++){ - fscanf(ficpar," %lf",¶m[i][j][k]); - if(mle==1){ - printf(" %lf",param[i][j][k]); - fprintf(ficlog," %lf",param[i][j][k]); - } - else - fprintf(ficlog," %lf",param[i][j][k]); - fprintf(ficparo," %lf",param[i][j][k]); - } - fscanf(ficpar,"\n"); + else{ + /* Read guess parameters */ + /* Reads comments: lines beginning with '#' */ + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); numlinepar++; - if(mle==1) - printf("\n"); - fprintf(ficlog,"\n"); - fprintf(ficparo,"\n"); + puts(line); + fputs(line,ficparo); + fputs(line,ficlog); } - } - fflush(ficlog); - - npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ + 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); - p=param[1][1]; - - /* Reads comments: lines beginning with '#' */ - while((c=getc(ficpar))=='#' && c!= EOF){ + p=param[1][1]; + + /* Reads comments: lines beginning with '#' */ + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + numlinepar++; + puts(line); + fputs(line,ficparo); + fputs(line,ficlog); + } ungetc(c,ficpar); - fgets(line, MAXLINE, ficpar); - numlinepar++; - puts(line); - fputs(line,ficparo); - fputs(line,ficlog); - } - ungetc(c,ficpar); - delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); - /* delti=vector(1,npar); *//* Scale of each paramater (output from hesscov) */ - for(i=1; i <=nlstate; i++){ - for(j=1; j <=nlstate+ndeath-1; j++){ - fscanf(ficpar,"%1d%1d",&i1,&j1); - if ((i1-i)*(j1-j)!=0){ - printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); - exit(1); - } - printf("%1d%1d",i,j); - fprintf(ficparo,"%1d%1d",i1,j1); - fprintf(ficlog,"%1d%1d",i1,j1); - for(k=1; k<=ncovmodel;k++){ - fscanf(ficpar,"%le",&delti3[i][j][k]); - printf(" %le",delti3[i][j][k]); - fprintf(ficparo," %le",delti3[i][j][k]); - fprintf(ficlog," %le",delti3[i][j][k]); + for(i=1; i <=nlstate; i++){ + for(j=1; j <=nlstate+ndeath-1; j++){ + fscanf(ficpar,"%1d%1d",&i1,&j1); + if ((i1-i)*(j1-j)!=0){ + printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); + exit(1); + } + printf("%1d%1d",i,j); + fprintf(ficparo,"%1d%1d",i1,j1); + fprintf(ficlog,"%1d%1d",i1,j1); + for(k=1; k<=ncovmodel;k++){ + fscanf(ficpar,"%le",&delti3[i][j][k]); + printf(" %le",delti3[i][j][k]); + fprintf(ficparo," %le",delti3[i][j][k]); + fprintf(ficlog," %le",delti3[i][j][k]); + } + fscanf(ficpar,"\n"); + numlinepar++; + printf("\n"); + fprintf(ficparo,"\n"); + fprintf(ficlog,"\n"); } - fscanf(ficpar,"\n"); - numlinepar++; - printf("\n"); - fprintf(ficparo,"\n"); - fprintf(ficlog,"\n"); } - } - fflush(ficlog); + fflush(ficlog); - delti=delti3[1][1]; + delti=delti3[1][1]; - /* free_ma3x(delti3,1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */ /* Hasn't to to freed here otherwise delti is no more allocated */ + /* free_ma3x(delti3,1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */ /* Hasn't to to freed here otherwise delti is no more allocated */ - /* Reads comments: lines beginning with '#' */ - while((c=getc(ficpar))=='#' && c!= EOF){ + /* Reads comments: lines beginning with '#' */ + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + numlinepar++; + puts(line); + fputs(line,ficparo); + fputs(line,ficlog); + } ungetc(c,ficpar); - fgets(line, MAXLINE, ficpar); - numlinepar++; - puts(line); - fputs(line,ficparo); - fputs(line,ficlog); - } - ungetc(c,ficpar); - matcov=matrix(1,npar,1,npar); - for(i=1; i <=npar; i++){ - fscanf(ficpar,"%s",&str); - if(mle==1) - printf("%s",str); - fprintf(ficlog,"%s",str); - fprintf(ficparo,"%s",str); - for(j=1; j <=i; j++){ - fscanf(ficpar," %le",&matcov[i][j]); - if(mle==1){ - printf(" %.5le",matcov[i][j]); + matcov=matrix(1,npar,1,npar); + for(i=1; i <=npar; i++){ + fscanf(ficpar,"%s",&str); + if(mle==1) + printf("%s",str); + fprintf(ficlog,"%s",str); + fprintf(ficparo,"%s",str); + for(j=1; j <=i; j++){ + fscanf(ficpar," %le",&matcov[i][j]); + if(mle==1){ + printf(" %.5le",matcov[i][j]); + } + fprintf(ficlog," %.5le",matcov[i][j]); + fprintf(ficparo," %.5le",matcov[i][j]); } - fprintf(ficlog," %.5le",matcov[i][j]); - fprintf(ficparo," %.5le",matcov[i][j]); + fscanf(ficpar,"\n"); + numlinepar++; + if(mle==1) + printf("\n"); + fprintf(ficlog,"\n"); + fprintf(ficparo,"\n"); } - fscanf(ficpar,"\n"); - numlinepar++; + for(i=1; i <=npar; i++) + for(j=i+1;j<=npar;j++) + matcov[i][j]=matcov[j][i]; + if(mle==1) printf("\n"); fprintf(ficlog,"\n"); - fprintf(ficparo,"\n"); - } - for(i=1; i <=npar; i++) - for(j=i+1;j<=npar;j++) - matcov[i][j]=matcov[j][i]; - - if(mle==1) - printf("\n"); - fprintf(ficlog,"\n"); - - fflush(ficlog); - - /*-------- Rewriting paramater file ----------*/ - strcpy(rfileres,"r"); /* "Rparameterfile */ - strcat(rfileres,optionfilefiname); /* Parameter file first name*/ - strcat(rfileres,"."); /* */ - strcat(rfileres,optionfilext); /* Other files have txt extension */ - if((ficres =fopen(rfileres,"w"))==NULL) { - printf("Problem writing new parameter file: %s\n", fileres);goto end; - fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end; - } - fprintf(ficres,"#%s\n",version); + fflush(ficlog); + + /*-------- Rewriting parameter file ----------*/ + strcpy(rfileres,"r"); /* "Rparameterfile */ + strcat(rfileres,optionfilefiname); /* Parameter file first name*/ + strcat(rfileres,"."); /* */ + strcat(rfileres,optionfilext); /* Other files have txt extension */ + if((ficres =fopen(rfileres,"w"))==NULL) { + printf("Problem writing new parameter file: %s\n", fileres);goto end; + fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end; + } + fprintf(ficres,"#%s\n",version); + } /* End of mle != -3 */ + /*-------- data file ----------*/ if((fic=fopen(datafile,"r"))==NULL) { - printf("Problem with datafile: %s\n", datafile);goto end; - fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end; + printf("Problem while opening datafile: %s\n", datafile);goto end; + fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);goto end; } n= lastobs; @@ -4252,40 +4727,117 @@ int main(int argc, char *argv[]) ncodemax=ivector(1,8); i=1; - while (fgets(line, MAXLINE, fic) != NULL) { - if ((i >= firstobs) && (i <=lastobs)) { - - for (j=maxwav;j>=1;j--){ - cutv(stra, strb,line,' '); s[j][i]=atoi(strb); - strcpy(line,stra); - cutv(stra, strb,line,'/'); anint[j][i]=(double)(atoi(strb)); strcpy(line,stra); - cutv(stra, strb,line,' '); mint[j][i]=(double)(atoi(strb)); strcpy(line,stra); + linei=0; + while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) { + linei=linei+1; + for(j=strlen(line); j>=0;j--){ /* Untabifies line */ + if(line[j] == '\t') + line[j] = ' '; + } + for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10)||(line[j]==13);j--){ + ; + }; + line[j+1]=0; /* Trims blanks at end of line */ + if(line[0]=='#'){ + fprintf(ficlog,"Comment line\n%s\n",line); + printf("Comment line\n%s\n",line); + continue; + } + + for (j=maxwav;j>=1;j--){ + cutv(stra, strb,line,' '); + errno=0; + lval=strtol(strb,&endptr,10); + /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/ + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%d' at line number %d %s for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav); + exit(1); } - - cutv(stra, strb,line,'/'); andc[i]=(double)(atoi(strb)); strcpy(line,stra); - cutv(stra, strb,line,' '); moisdc[i]=(double)(atoi(strb)); strcpy(line,stra); - - cutv(stra, strb,line,'/'); annais[i]=(double)(atoi(strb)); strcpy(line,stra); - cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra); - - cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra); - for (j=ncovcol;j>=1;j--){ - cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra); - } - lstra=strlen(stra); - if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */ - stratrunc = &(stra[lstra-9]); - num[i]=atol(stratrunc); + s[j][i]=lval; + + strcpy(line,stra); + cutv(stra, strb,line,' '); + if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ + } + else if(iout=sscanf(strb,"%s.") != 0){ + month=99; + year=9999; + }else{ + printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j); + exit(1); } - else - num[i]=atol(stra); - - /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){ - printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/ - - i=i+1; + anint[j][i]= (double) year; + mint[j][i]= (double)month; + strcpy(line,stra); + } /* ENd Waves */ + + cutv(stra, strb,line,' '); + if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ } - } + else if(iout=sscanf(strb,"%s.",dummy) != 0){ + month=99; + year=9999; + }else{ + printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line); + exit(1); + } + andc[i]=(double) year; + moisdc[i]=(double) month; + strcpy(line,stra); + + cutv(stra, strb,line,' '); + if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ + } + else if(iout=sscanf(strb,"%s.") != 0){ + month=99; + year=9999; + }else{ + printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j); + exit(1); + } + annais[i]=(double)(year); + moisnais[i]=(double)(month); + strcpy(line,stra); + + cutv(stra, strb,line,' '); + errno=0; + lval=strtol(strb,&endptr,10); + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a weight. Exiting.\n",lval, i,line,linei); + exit(1); + } + weight[i]=(double)(lval); + strcpy(line,stra); + + for (j=ncovcol;j>=1;j--){ + cutv(stra, strb,line,' '); + errno=0; + lval=strtol(strb,&endptr,10); + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a covar (meaning 0 for the reference or 1). Exiting.\n",lval, linei,i, line); + exit(1); + } + if(lval <-1 || lval >1){ + printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a value of the %d covar (meaning 0 for the reference or 1. IMaCh does not build design variables, do it your self). Exiting.\n",lval,linei, i,line,j); + exit(1); + } + covar[j][i]=(double)(lval); + strcpy(line,stra); + } + lstra=strlen(stra); + + if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */ + stratrunc = &(stra[lstra-9]); + num[i]=atol(stratrunc); + } + else + num[i]=atol(stra); + /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){ + printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/ + + i=i+1; + } /* End loop reading data */ + fclose(fic); /* printf("ii=%d", ij); scanf("%d",i);*/ imx=i-1; /* Number of individuals */ @@ -4299,12 +4851,12 @@ int main(int argc, char *argv[]) if (s[4][i]==9) s[4][i]=-1; printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i]));}*/ - for (i=1; i<=imx; i++) + /* for (i=1; i<=imx; i++) */ /*if ((s[3][i]==3) || (s[4][i]==3)) weight[i]=0.08; else weight[i]=1;*/ - /* Calculation of the number of parameter from char model*/ + /* Calculation of the number of parameters from char model */ Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */ Tprod=ivector(1,15); Tvaraff=ivector(1,15); @@ -4383,8 +4935,7 @@ int main(int argc, char *argv[]) printf("cptcovprod=%d ", cptcovprod); fprintf(ficlog,"cptcovprod=%d ", cptcovprod); - scanf("%d ",i); - fclose(fic);*/ + scanf("%d ",i);*/ /* if(mle==1){*/ if (weightopt != 1) { /* Maximisation without weights*/ @@ -4417,7 +4968,7 @@ int main(int argc, char *argv[]) for (i=1; i<=imx; i++) { agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); for(m=firstpass; (m<= lastpass); m++){ - if(s[m][i] >0){ + if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){ if (s[m][i] >= nlstate+1) { if(agedc[i]>0) if((int)moisdc[i]!=99 && (int)andc[i]!=9999) @@ -4433,8 +4984,7 @@ int main(int argc, char *argv[]) } } else if(s[m][i] !=9){ /* Standard case, age in fractional - years but with the precision of a - month */ + years but with the precision of a month */ agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]); if((int)mint[m][i]==99 || (int)anint[m][i]==9999) agev[m][i]=1; @@ -4477,9 +5027,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); @@ -4533,7 +5085,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); } @@ -4546,6 +5101,8 @@ int main(int argc, char *argv[]) /*--------- index.htm --------*/ strcpy(optionfilehtm,optionfilefiname); /* Main html file */ + if(mle==-3) + strcat(optionfilehtm,"-mort"); strcat(optionfilehtm,".htm"); if((fichtm=fopen(optionfilehtm,"w"))==NULL) { printf("Problem with %s \n",optionfilehtm), exit(0); @@ -4569,11 +5126,13 @@ Title=%s
                Datafile=%s Firstpass=%d La \n\
                \

                • Parameter files

                  \n\ + - Parameter file: %s.%s
                  \n\ - Copy of the parameter file: o%s
                  \n\ - 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,\ + optionfilefiname,optionfilext,optionfilefiname,optionfilext,\ fileres,fileres,\ filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart); fflush(fichtm); @@ -4581,13 +5140,10 @@ Title=%s
                Datafile=%s Firstpass=%d La 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'. */ - freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart); fprintf(fichtm,"\n"); fprintf(fichtm,"
                Total number of observations=%d
                \n\ @@ -4606,33 +5162,154 @@ 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 (m=firstpass; m<=lastpass; m++) + if (s[m][i]>nlstate) { + dcwave[i]=m; + /* 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; + }else 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; + } - /*--------- 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); - + for (i=1;i<=NDIM;i++) { + for (j=1;j<=NDIM;j++) + ximort[i][j]=(i == j ? 1.0 : 0.0); + } + + p[1]=0.0268; p[NDIM]=0.083; + /*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); + + 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])); + + lsurv=vector(1,AGESUP); + lpop=vector(1,AGESUP); + tpop=vector(1,AGESUP); + lsurv[agegomp]=100000; + + for (k=agegomp;k<=AGESUP;k++) { + agemortsup=k; + if (p[1]*exp(p[2]*(k-agegomp))>1) break; + } + + for (k=agegomp;k=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); @@ -4646,473 +5323,559 @@ 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,pathcd); /* 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 (period or stable prevalence) --------------*/ + + strcpy(filerespl,"pl"); + strcat(filerespl,fileres); + if((ficrespl=fopen(filerespl,"w"))==NULL) { + printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end; + fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end; + } + printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl); + fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl); + pstamp(ficrespl); + fprintf(ficrespl,"# Period (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 */ - - /* hstepm=1; aff par mois*/ + agelim=AGESUP; + hstepm=stepsize*YEARM; /* Every year of age */ + hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ - 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"); + /* hstepm=1; aff par mois*/ + pstamp(ficrespij); + 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,strstart); - 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 ------------*/ - - strcpy(filerest,"t"); - strcat(filerest,fileres); - if((ficrest=fopen(filerest,"w"))==NULL) { - printf("Problem with total LE resultfile: %s\n", filerest);goto end; - fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end; - } - printf("Computing Total LEs with variances: file '%s' \n", filerest); - fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest); - - - strcpy(filerese,"e"); - strcat(filerese,fileres); - if((ficreseij=fopen(filerese,"w"))==NULL) { - printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0); - fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0); - } - printf("Computing Health Expectancies: result on file '%s' \n", filerese); - fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese); - - strcpy(fileresv,"v"); - strcat(fileresv,fileres); - if((ficresvij=fopen(fileresv,"w"))==NULL) { - printf("Problem with variance resultfile: %s\n", fileresv);exit(0); - fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0); - } - printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); - fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); - - /* 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); - */ + /*---------- 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 Life expectancies with their standard errors: file '%s' \n", filerest); + fprintf(ficlog,"Computing Total Life expectancies with their standard errors: 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(fileresstde,"stde"); + strcat(fileresstde,fileres); + if((ficresstdeij=fopen(fileresstde,"w"))==NULL) { + printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0); + fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0); + } + printf("Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); + fprintf(ficlog,"Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); + + strcpy(filerescve,"cve"); + strcat(filerescve,fileres); + if((ficrescveij=fopen(filerescve,"w"))==NULL) { + printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0); + fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0); + } + printf("Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); + fprintf(ficlog,"Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); + + 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); + */ - 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); - - 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); - } + 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#****** "); + fprintf(ficresstdeij,"\n#****** "); + fprintf(ficrescveij,"\n#****** "); + for(j=1;j<=cptcoveff;j++) { + fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + } + fprintf(ficreseij,"******\n"); + fprintf(ficresstdeij,"******\n"); + fprintf(ficrescveij,"******\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, strstart); + cvevsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart); - 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]);*/ - } - epj[nlstate+1] +=epj[j]; + 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, strstart); + 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, strstart); } - 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])); - } + pstamp(ficrest); + fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# Age ( e.. (std) "); + for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); 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"); + 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]);*/ + } + 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])); + } + 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(ficresstdeij); + fclose(ficrescveij); + fclose(ficresvij); + fclose(ficrest); + fclose(ficpar); + + /*------- Variance of period (stable) prevalence------*/ + + strcpy(fileresvpl,"vpl"); + strcat(fileresvpl,fileres); + if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { + printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl); + exit(0); + } + printf("Computing Variance-covariance of period (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,strstart); + 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); - 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); + /*---------- 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(prlim,1,nlstate,1,nlstate); + 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_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX); + free_imatrix(codtab,1,100,1,10); fflush(fichtm); fflush(ficgp); if((nberr >0) || (nbwarn>0)){ - printf("End of Imach with %d errors and/or warnings %d\n",nberr,nbwarn); + 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"); @@ -5124,8 +5887,8 @@ ageminpar, agemax, s[lastpass][imx], age 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\n",strstart, strtend); + printf("Local time at start %s\nLocal time 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); @@ -5140,13 +5903,38 @@ ageminpar, agemax, s[lastpass][imx], age fclose(ficlog); /*------ End -----------*/ - chdir(path); - strcpy(plotcmd,GNUPLOTPROGRAM); - strcat(plotcmd," "); - strcat(plotcmd,optionfilegnuplot); - printf("Starting graphs with: %s",plotcmd);fflush(stdout); + + printf("Before Current directory %s!\n",pathcd); + if(chdir(pathcd) != 0) + printf("Can't move to directory %s!\n",path); + if(getcwd(pathcd,MAXLINE) > 0) + printf("Current directory %s!\n",pathcd); + /*strcat(plotcmd,CHARSEPARATOR);*/ + sprintf(plotcmd,"gnuplot"); +#ifndef UNIX + sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach); +#endif + if(!stat(plotcmd,&info)){ + printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout); + if(!stat(getenv("GNUPLOTBIN"),&info)){ + printf("Error gnuplot program not found: %s Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout); + }else + strcpy(pplotcmd,plotcmd); +#ifdef UNIX + strcpy(plotcmd,GNUPLOTPROGRAM); + if(!stat(plotcmd,&info)){ + printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout); + }else + strcpy(pplotcmd,plotcmd); +#endif + }else + strcpy(pplotcmd,plotcmd); + + sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot); + printf("Starting graphs with: %s\n",plotcmd);fflush(stdout); + if((outcmd=system(plotcmd)) != 0){ - printf(" Problem with gnuplot\n"); + printf("\n Problem with gnuplot\n"); } printf(" Wait..."); while (z[0] != 'q') { @@ -5154,7 +5942,10 @@ ageminpar, agemax, s[lastpass][imx], age printf("\nType e to edit output files, g to graph again and q for exiting: "); scanf("%s",z); /* if (z[0] == 'c') system("./imach"); */ - if (z[0] == 'e') system(optionfilehtm); + 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); }