|
|
| version 1.99, 2004/06/05 08:57:40 | version 1.117, 2006/03/14 17:16:22 |
|---|---|
| Line 1 | Line 1 |
| /* $Id$ | /* $Id$ |
| $State$ | $State$ |
| $Log$ | $Log$ |
| 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 | Revision 1.99 2004/06/05 08:57:40 brouard |
| *** empty log message *** | *** empty log message *** |
| Line 146 | Line 217 |
| hPijx. | hPijx. |
| Also this programme outputs the covariance matrix of the parameters but also | 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). | Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). |
| Institut national d'études démographiques, Paris. | Institut national d'études démographiques, Paris. |
| Line 175 | Line 246 |
| begin-prev-date,... | begin-prev-date,... |
| open gnuplot file | open gnuplot file |
| open html file | open html file |
| stable prevalence | period (stable) prevalence |
| for age prevalim() | for age prevalim() |
| h Pij x | h Pij x |
| variance of p varprob | variance of p varprob |
| Line 187 | Line 258 |
| varevsij() | varevsij() |
| if popbased==1 varevsij(,popbased) | if popbased==1 varevsij(,popbased) |
| total life expectancies | total life expectancies |
| Variance of stable prevalence | Variance of period (stable) prevalence |
| end | end |
| */ | */ |
| Line 197 | Line 268 |
| #include <math.h> | #include <math.h> |
| #include <stdio.h> | #include <stdio.h> |
| #include <stdlib.h> | #include <stdlib.h> |
| #include <string.h> | |
| #include <unistd.h> | #include <unistd.h> |
| #include <limits.h> | |
| #include <sys/types.h> | |
| #include <sys/stat.h> | |
| #include <errno.h> | |
| extern int errno; | |
| /* #include <sys/time.h> */ | /* #include <sys/time.h> */ |
| #include <time.h> | #include <time.h> |
| #include "timeval.h" | #include "timeval.h" |
| Line 207 | Line 285 |
| /* #define _(String) gettext (String) */ | /* #define _(String) gettext (String) */ |
| #define MAXLINE 256 | #define MAXLINE 256 |
| #define GNUPLOTPROGRAM "gnuplot" | #define GNUPLOTPROGRAM "gnuplot" |
| /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ | /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ |
| #define FILENAMELENGTH 132 | #define FILENAMELENGTH 132 |
| /*#define DEBUG*/ | |
| /*#define windows*/ | |
| #define GLOCK_ERROR_NOPATH -1 /* empty path */ | #define GLOCK_ERROR_NOPATH -1 /* empty path */ |
| #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ | #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ |
| Line 227 | Line 305 |
| #define AGESUP 130 | #define AGESUP 130 |
| #define AGEBASE 40 | #define AGEBASE 40 |
| #define AGEGOMP 10. /* Minimal age for Gompertz adjustment */ | #define AGEGOMP 10. /* Minimal age for Gompertz adjustment */ |
| #ifdef unix | #ifdef UNIX |
| #define DIRSEPARATOR '/' | #define DIRSEPARATOR '/' |
| #define CHARSEPARATOR "/" | |
| #define ODIRSEPARATOR '\\' | #define ODIRSEPARATOR '\\' |
| #else | #else |
| #define DIRSEPARATOR '\\' | #define DIRSEPARATOR '\\' |
| #define CHARSEPARATOR "\\" | |
| #define ODIRSEPARATOR '/' | #define ODIRSEPARATOR '/' |
| #endif | #endif |
| /* $Id$ */ | /* $Id$ */ |
| /* $State$ */ | /* $State$ */ |
| char version[]="Imach version 0.97b, May 2004, INED-EUROREVES "; | char version[]="Imach version 0.98d, March 2006, INED-EUROREVES-Institut de longevite "; |
| char fullversion[]="$Revision$ $Date$"; | char fullversion[]="$Revision$ $Date$"; |
| char strstart[80]; | |
| char optionfilext[10], optionfilefiname[FILENAMELENGTH]; | |
| int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ | int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ |
| int nvar; | int nvar; |
| int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; | int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; |
| Line 252 int popbased=0; | Line 334 int popbased=0; |
| int *wav; /* Number of waves for this individuual 0 is possible */ | int *wav; /* Number of waves for this individuual 0 is possible */ |
| int maxwav; /* Maxim number of waves */ | int maxwav; /* Maxim number of waves */ |
| int jmin, jmax; /* min, max spacing between 2 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 | int gipmx, gsw; /* Global variables on the number of contributions |
| to the likelihood and the sum of weights (done by funcone)*/ | to the likelihood and the sum of weights (done by funcone)*/ |
| int mle, weightopt; | int mle, weightopt; |
| Line 262 int **bh; /* bh[mi][i] is the bias (+ or | Line 345 int **bh; /* bh[mi][i] is the bias (+ or |
| double jmean; /* Mean space between 2 waves */ | double jmean; /* Mean space between 2 waves */ |
| double **oldm, **newm, **savm; /* Working pointers to matrices */ | double **oldm, **newm, **savm; /* Working pointers to matrices */ |
| double **oldms, **newms, **savms; /* Fixed 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; | FILE *ficlog, *ficrespow; |
| int globpr; /* Global variable for printing or not */ | int globpr; /* Global variable for printing or not */ |
| double fretone; /* Only one call to likelihood */ | double fretone; /* Only one call to likelihood */ |
| Line 276 FILE *ficresprobmorprev; | Line 359 FILE *ficresprobmorprev; |
| FILE *fichtm, *fichtmcov; /* Html File */ | FILE *fichtm, *fichtmcov; /* Html File */ |
| FILE *ficreseij; | FILE *ficreseij; |
| char filerese[FILENAMELENGTH]; | char filerese[FILENAMELENGTH]; |
| FILE *ficresstdeij; | |
| char fileresstde[FILENAMELENGTH]; | |
| FILE *ficrescveij; | |
| char filerescve[FILENAMELENGTH]; | |
| FILE *ficresvij; | FILE *ficresvij; |
| char fileresv[FILENAMELENGTH]; | char fileresv[FILENAMELENGTH]; |
| FILE *ficresvpl; | FILE *ficresvpl; |
| char fileresvpl[FILENAMELENGTH]; | char fileresvpl[FILENAMELENGTH]; |
| char title[MAXLINE]; | char title[MAXLINE]; |
| char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH]; | char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH]; |
| char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH]; | char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH]; |
| char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH]; | char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH]; |
| char command[FILENAMELENGTH]; | char command[FILENAMELENGTH]; |
| int outcmd=0; | int outcmd=0; |
| Line 304 long time_value; | Line 391 long time_value; |
| extern long time(); | extern long time(); |
| char strcurr[80], strfor[80]; | char strcurr[80], strfor[80]; |
| char *endptr; | |
| long lval; | |
| #define NR_END 1 | #define NR_END 1 |
| #define FREE_ARG char* | #define FREE_ARG char* |
| #define FTOL 1.0e-10 | #define FTOL 1.0e-10 |
| Line 352 double *weight; | Line 442 double *weight; |
| int **s; /* Status */ | int **s; /* Status */ |
| double *agedc, **covar, idx; | double *agedc, **covar, idx; |
| int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff; | int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff; |
| double *lsurv, *lpop, *tpop; | |
| double ftol=FTOL; /* Tolerance for computing Max Likelihood */ | double ftol=FTOL; /* Tolerance for computing Max Likelihood */ |
| double ftolhess; /* Tolerance for computing hessian */ | double ftolhess; /* Tolerance for computing hessian */ |
| Line 359 double ftolhess; /* Tolerance for comput | Line 450 double ftolhess; /* Tolerance for comput |
| /**************** split *************************/ | /**************** split *************************/ |
| static int split( char *path, char *dirc, char *name, char *ext, char *finame ) | 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) | /* 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) | the name of the file (name), its extension only (ext) and its first part of the name (finame) |
| */ | */ |
| char *ss; /* pointer */ | char *ss; /* pointer */ |
| Line 368 static int split( char *path, char *dirc | Line 459 static int split( char *path, char *dirc |
| l1 = strlen(path ); /* length of path */ | l1 = strlen(path ); /* length of path */ |
| if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH ); | if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH ); |
| ss= strrchr( path, DIRSEPARATOR ); /* find last / */ | 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) | /*if(strrchr(path, ODIRSEPARATOR )==NULL) |
| printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/ | printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/ |
| /* get current working directory */ | /* get current working directory */ |
| Line 376 static int split( char *path, char *dirc | Line 468 static int split( char *path, char *dirc |
| if ( getcwd( dirc, FILENAME_MAX ) == NULL ) { | if ( getcwd( dirc, FILENAME_MAX ) == NULL ) { |
| return( GLOCK_ERROR_GETCWD ); | 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 */ | } else { /* strip direcotry from path */ |
| ss++; /* after this, the filename */ | ss++; /* after this, the filename */ |
| l2 = strlen( ss ); /* length of filename */ | l2 = strlen( ss ); /* length of filename */ |
| Line 384 static int split( char *path, char *dirc | Line 477 static int split( char *path, char *dirc |
| strcpy( name, ss ); /* save file name */ | strcpy( name, ss ); /* save file name */ |
| strncpy( dirc, path, l1 - l2 ); /* now the directory */ | strncpy( dirc, path, l1 - l2 ); /* now the directory */ |
| dirc[l1-l2] = 0; /* add zero */ | 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 */ | l1 = strlen( dirc ); /* length of directory */ |
| /*#ifdef windows | if( dirc[l1-1] != DIRSEPARATOR ){ |
| if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; } | dirc[l1] = DIRSEPARATOR; |
| #else | dirc[l1+1] = 0; |
| if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; } | printf(" DIRC3 = %s \n",dirc); |
| #endif | } |
| */ | |
| ss = strrchr( name, '.' ); /* find last / */ | ss = strrchr( name, '.' ); /* find last / */ |
| if (ss >0){ | if (ss >0){ |
| ss++; | ss++; |
| Line 401 static int split( char *path, char *dirc | Line 495 static int split( char *path, char *dirc |
| strncpy( finame, name, l1-l2); | strncpy( finame, name, l1-l2); |
| finame[l1-l2]= 0; | finame[l1-l2]= 0; |
| } | } |
| return( 0 ); /* we're done */ | return( 0 ); /* we're done */ |
| } | } |
| Line 433 int nbocc(char *s, char occ) | Line 528 int nbocc(char *s, char occ) |
| void cutv(char *u,char *v, char*t, 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 | /* cuts string t into u and v where u ends before first occurence of char 'occ' |
| and v is after occ excluding it too : ex cutv(u,v,"abcdef2ghi2j",2) | and v starts after first occurence of char 'occ' : ex cutv(u,v,"abcdef2ghi2j",'2') |
| gives u="abcedf" and v="ghi2j" */ | gives u="abcedf" and v="ghi2j" */ |
| int i,lg,j,p=0; | int i,lg,j,p=0; |
| i=0; | i=0; |
| Line 876 void powell(double p[], double **xi, int | Line 971 void powell(double p[], double **xi, int |
| fprintf(ficrespow,"\n");fflush(ficrespow); | fprintf(ficrespow,"\n");fflush(ficrespow); |
| if(*iter <=3){ | if(*iter <=3){ |
| tm = *localtime(&curr_time.tv_sec); | tm = *localtime(&curr_time.tv_sec); |
| strcpy(strcurr,asctime(&tmf)); | strcpy(strcurr,asctime(&tm)); |
| /* asctime_r(&tm,strcurr); */ | /* asctime_r(&tm,strcurr); */ |
| forecast_time=curr_time; | forecast_time=curr_time; |
| itmp = strlen(strcurr); | itmp = strlen(strcurr); |
| if(strcurr[itmp-1]=='\n') | if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */ |
| strcurr[itmp-1]='\0'; | 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); | 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); | fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec); |
| Line 892 void powell(double p[], double **xi, int | Line 987 void powell(double p[], double **xi, int |
| itmp = strlen(strfor); | itmp = strlen(strfor); |
| if(strfor[itmp-1]=='\n') | if(strfor[itmp-1]=='\n') |
| strfor[itmp-1]='\0'; | 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); | 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 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 i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); |
| } | } |
| } | } |
| for (i=1;i<=n;i++) { | for (i=1;i<=n;i++) { |
| Line 987 void powell(double p[], double **xi, int | Line 1082 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) | double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij) |
| { | { |
| Line 1232 double func( double *x) | Line 1327 double func( double *x) |
| } /* end mult */ | } /* end mult */ |
| /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ | /*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. | /* But now since version 0.9 we anticipate for bias at large stepm. |
| * If stepm is larger than one month (smallest stepm) and if the exact delay | * 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 | * (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 | * 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 | * 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 | * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the |
| * probability in order to take into account the bias as a fraction of the way | * 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 | * from savm to out if bh is negative or even beyond if bh is positive. bh varies |
| * -stepm/2 to stepm/2 . | * -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 the same as for previous versions of Imach. |
| * For stepm > 1 the results are less biased than in previous versions. | * For stepm > 1 the results are less biased than in previous versions. |
| Line 1247 double func( double *x) | Line 1342 double func( double *x) |
| s1=s[mw[mi][i]][i]; | s1=s[mw[mi][i]][i]; |
| s2=s[mw[mi+1][i]][i]; | s2=s[mw[mi+1][i]][i]; |
| bbh=(double)bh[mi][i]/(double)stepm; | 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. | * 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]));*/ | /* 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){ | if( s2 > nlstate){ |
| /* i.e. if s2 is a death state and if the date of death is known then the contribution | /* i.e. if s2 is a death state and if the date of death is known |
| to the likelihood is the probability to die between last step unit time and current | then the contribution to the likelihood is the probability to |
| step unit time, which is also the differences between probability to die before dh | die between last step unit time and current step unit time, |
| and probability to die before dh-stepm . | 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 | In version up to 0.92 likelihood was computed |
| as if date of death was unknown. Death was treated as any other | as if date of death was unknown. Death was treated as any other |
| health state: the date of the interview describes the actual state | health state: the date of the interview describes the actual state |
| Line 1275 double func( double *x) | Line 1371 double func( double *x) |
| lower mortality. | lower mortality. |
| */ | */ |
| lli=log(out[s1][s2] - savm[s1][s2]); | 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= 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= 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 */ | /* 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 */ |
| } | } |
| Line 1571 void mlikeli(FILE *ficres,double p[], in | Line 1688 void mlikeli(FILE *ficres,double p[], in |
| powell(p,xi,npar,ftol,&iter,&fret,func); | powell(p,xi,npar,ftol,&iter,&fret,func); |
| free_matrix(xi,1,npar,1,npar); | |
| fclose(ficrespow); | fclose(ficrespow); |
| printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p)); | 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)); | fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); |
| Line 1838 void lubksb(double **a, int n, int *indx | Line 1956 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 ********************/ | /************ 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 */ | { /* Some frequencies */ |
| int i, m, jk, k1,i1, j1, bool, z1,z2,j; | int i, m, jk, k1,i1, j1, bool, z1,z2,j; |
| Line 1847 void freqsummary(char fileres[], int ia | Line 1970 void freqsummary(char fileres[], int ia |
| double ***freq; /* Frequencies */ | double ***freq; /* Frequencies */ |
| double *pp, **prop; | double *pp, **prop; |
| double pos,posprop, k2, dateintsum=0,k2cpt=0; | double pos,posprop, k2, dateintsum=0,k2cpt=0; |
| FILE *ficresp; | |
| char fileresp[FILENAMELENGTH]; | char fileresp[FILENAMELENGTH]; |
| pp=vector(1,nlstate); | pp=vector(1,nlstate); |
| Line 1859 void freqsummary(char fileres[], int ia | Line 1981 void freqsummary(char fileres[], int ia |
| fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); | fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); |
| exit(0); | 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; | j1=0; |
| j=cptcoveff; | j=cptcoveff; |
| Line 1872 void freqsummary(char fileres[], int ia | Line 1994 void freqsummary(char fileres[], int ia |
| j1++; | j1++; |
| /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); | /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
| scanf("%d", i);*/ | scanf("%d", i);*/ |
| for (i=-1; i<=nlstate+ndeath; i++) | for (i=-5; i<=nlstate+ndeath; i++) |
| for (jk=-1; jk<=nlstate+ndeath; jk++) | for (jk=-5; jk<=nlstate+ndeath; jk++) |
| for(m=iagemin; m <= iagemax+3; m++) | for(m=iagemin; m <= iagemax+3; m++) |
| freq[i][jk][m]=0; | freq[i][jk][m]=0; |
| Line 1912 void freqsummary(char fileres[], int ia | Line 2034 void freqsummary(char fileres[], int ia |
| } | } |
| /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ | /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
| pstamp(ficresp); | |
| if (cptcovn>0) { | if (cptcovn>0) { |
| fprintf(ficresp, "\n#********** Variable "); | fprintf(ficresp, "\n#********** Variable "); |
| for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); | for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); |
| Line 1998 void freqsummary(char fileres[], int ia | Line 2120 void freqsummary(char fileres[], int ia |
| dateintmean=dateintsum/k2cpt; | dateintmean=dateintsum/k2cpt; |
| fclose(ficresp); | 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_vector(pp,1,nlstate); |
| free_matrix(prop,1,nlstate,iagemin, iagemax+3); | free_matrix(prop,1,nlstate,iagemin, iagemax+3); |
| /* End of Freq */ | /* End of Freq */ |
| Line 2107 void concatwav(int wav[], int **dh, int | Line 2229 void concatwav(int wav[], int **dh, int |
| mi=0; | mi=0; |
| m=firstpass; | m=firstpass; |
| while(s[m][i] <= nlstate){ | 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; | mw[++mi][i]=m; |
| if(m >=lastpass) | if(m >=lastpass) |
| break; | break; |
| Line 2125 void concatwav(int wav[], int **dh, int | Line 2247 void concatwav(int wav[], int **dh, int |
| if(mi==0){ | if(mi==0){ |
| nbwarn++; | nbwarn++; |
| if(first==0){ | 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; | first=1; |
| } | } |
| if(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 mi==0 */ |
| } /* End individuals */ | } /* End individuals */ |
| Line 2147 void concatwav(int wav[], int **dh, int | Line 2269 void concatwav(int wav[], int **dh, int |
| nberr++; | 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]); | 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 */ | 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,"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; | k=k+1; |
| if (j >= jmax) jmax=j; | if (j >= jmax){ |
| if (j <= jmin) jmin=j; | jmax=j; |
| ijmax=i; | |
| } | |
| if (j <= jmin){ | |
| jmin=j; | |
| ijmin=i; | |
| } | |
| sum=sum+j; | sum=sum+j; |
| /*if (j<0) printf("j=%d num=%d \n",j,i);*/ | /*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);*/ | /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ |
| Line 2161 void concatwav(int wav[], int **dh, int | Line 2289 void concatwav(int wav[], int **dh, int |
| } | } |
| else{ | else{ |
| j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); | 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; | k=k+1; |
| if (j >= jmax) jmax=j; | if (j >= jmax) { |
| else if (j <= jmin)jmin=j; | 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); */ | /* 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]);*/ | /*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){ | if(j<0){ |
| Line 2207 void concatwav(int wav[], int **dh, int | Line 2342 void concatwav(int wav[], int **dh, int |
| } /* end wave */ | } /* end wave */ |
| } | } |
| jmean=sum/k; | jmean=sum/k; |
| printf("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 Max=%d Mean=%f\n\n ",jmin, jmax,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 ****************************/ | /*********** Tricode ****************************/ |
| Line 2254 void tricode(int *Tvar, int **nbcode, in | Line 2389 void tricode(int *Tvar, int **nbcode, in |
| for (k=0; k< maxncov; k++) Ndum[k]=0; | for (k=0; k< maxncov; k++) Ndum[k]=0; |
| for (i=1; i<=ncovmodel-2; i++) { | 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]; | ij=Tvar[i]; |
| Ndum[ij]++; | Ndum[ij]++; |
| } | } |
| Line 2272 void tricode(int *Tvar, int **nbcode, in | Line 2407 void tricode(int *Tvar, int **nbcode, in |
| /*********** Health Expectancies ****************/ | /*********** 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. */ | |
| /* Computing Variances of health expectancies */ | |
| /* Gradient is computed with plus gp and minus gm. Code is duplicated in order to | |
| decrease memory allocation */ | |
| 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 */ | /* Covariances of health expectancies eij and of total life expectancies according |
| int i, j, nhstepm, hstepm, h, nstepm, k, cptj; | to initial status i, ei. . |
| */ | |
| int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji; | |
| double age, agelim, hf; | double age, agelim, hf; |
| double ***p3mat,***varhe; | double ***p3matp, ***p3matm, ***varhe; |
| double **dnewm,**doldm; | double **dnewm,**doldm; |
| double *xp; | double *xp, *xm; |
| double **gp, **gm; | double **gp, **gm; |
| double ***gradg, ***trgradg; | double ***gradg, ***trgradg; |
| int theta; | int theta; |
| double eip, vip; | |
| varhe=ma3x(1,nlstate*nlstate,1,nlstate*nlstate,(int) bage, (int) fage); | varhe=ma3x(1,nlstate*nlstate,1,nlstate*nlstate,(int) bage, (int) fage); |
| xp=vector(1,npar); | xp=vector(1,npar); |
| xm=vector(1,npar); | |
| dnewm=matrix(1,nlstate*nlstate,1,npar); | dnewm=matrix(1,nlstate*nlstate,1,npar); |
| doldm=matrix(1,nlstate*nlstate,1,nlstate*nlstate); | doldm=matrix(1,nlstate*nlstate,1,nlstate*nlstate); |
| fprintf(ficreseij,"# Health expectancies\n"); | pstamp(ficresstdeij); |
| fprintf(ficreseij,"# Age"); | fprintf(ficresstdeij,"# Health expectancies with standard errors\n"); |
| for(i=1; i<=nlstate;i++) | fprintf(ficresstdeij,"# Age"); |
| for(i=1; i<=nlstate;i++){ | |
| for(j=1; j<=nlstate;j++) | for(j=1; j<=nlstate;j++) |
| fprintf(ficreseij," %1d-%1d (SE)",i,j); | fprintf(ficresstdeij," e%1d%1d (SE)",i,j); |
| fprintf(ficreseij,"\n"); | 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){ | if(estepm < stepm){ |
| printf ("Problem %d lower than %d\n",estepm, stepm); | printf ("Problem %d lower than %d\n",estepm, stepm); |
| } | } |
| Line 2326 void evsij(char fileres[], double ***eij | Line 2583 void evsij(char fileres[], double ***eij |
| */ | */ |
| hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ | 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; | agelim=AGESUP; |
| for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ | nstepm=(int) rint((agelim-age)*YEARM/stepm); |
| /* nhstepm age range expressed in number of stepm */ | /* Typically if 20 years nstepm = 20*12/6=40 stepm */ |
| nstepm=(int) rint((agelim-age)*YEARM/stepm); | /* if (stepm >= YEARM) hstepm=1;*/ |
| /* Typically if 20 years nstepm = 20*12/6=40 stepm */ | nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ |
| /* 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); |
| p3mat=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); | gradg=ma3x(0,nhstepm,1,npar,1,nlstate*nlstate); |
| gp=matrix(0,nhstepm,1,nlstate*nlstate); | trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar); |
| gm=matrix(0,nhstepm,1,nlstate*nlstate); | 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 | /* Computed by stepm unit matrices, product of hstepm matrices, stored |
| in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ | 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. */ | hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ |
| /* Computing Variances of health expectancies */ | /* Computing Variances of health expectancies */ |
| /* Gradient is computed with plus gp and minus gm. Code is duplicated in order to | |
| for(theta=1; theta <=npar; theta++){ | decrease memory allocation */ |
| for(theta=1; theta <=npar; theta++){ | |
| for(i=1; i<=npar; i++){ | for(i=1; i<=npar; i++){ |
| xp[i] = x[i] + (i==theta ?delti[theta]:0); | 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(j=1; j<= nlstate; j++){ |
| for(i=1; i<=nlstate; i++){ | for(i=1; i<=nlstate; i++){ |
| cptj=cptj+1; | for(h=0; h<=nhstepm-1; h++){ |
| for(h=0, gp[h][cptj]=0.; h<=nhstepm-1; h++){ | gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.; |
| gp[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.; | gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.; |
| } | } |
| } | } |
| } | } |
| for(ij=1; ij<= nlstate*nlstate; ij++) | |
| 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(h=0; h<=nhstepm-1; h++){ | 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 */ |
| /* End theta */ | |
| for(h=0; h<=nhstepm-1; h++) | |
| trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar); | |
| for(h=0; h<=nhstepm-1; h++) | |
| for(j=1; j<=nlstate*nlstate;j++) | for(j=1; j<=nlstate*nlstate;j++) |
| for(theta=1; theta <=npar; theta++) | for(theta=1; theta <=npar; theta++) |
| trgradg[h][j][theta]=gradg[h][theta][j]; | trgradg[h][j][theta]=gradg[h][theta][j]; |
| for(i=1;i<=nlstate*nlstate;i++) | for(ij=1;ij<=nlstate*nlstate;ij++) |
| for(j=1;j<=nlstate*nlstate;j++) | for(ji=1;ji<=nlstate*nlstate;ji++) |
| varhe[i][j][(int)age] =0.; | varhe[ij][ji][(int)age] =0.; |
| printf("%d|",(int)age);fflush(stdout); | printf("%d|",(int)age);fflush(stdout); |
| fprintf(ficlog,"%d|",(int)age);fflush(ficlog); | fprintf(ficlog,"%d|",(int)age);fflush(ficlog); |
| Line 2404 void evsij(char fileres[], double ***eij | Line 2648 void evsij(char fileres[], double ***eij |
| for(k=0;k<=nhstepm-1;k++){ | for(k=0;k<=nhstepm-1;k++){ |
| matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov); | 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]); | matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]); |
| for(i=1;i<=nlstate*nlstate;i++) | for(ij=1;ij<=nlstate*nlstate;ij++) |
| for(j=1;j<=nlstate*nlstate;j++) | for(ji=1;ji<=nlstate*nlstate;ji++) |
| varhe[i][j][(int)age] += doldm[i][j]*hf*hf; | varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf; |
| } | } |
| } | } |
| /* Computing expectancies */ | /* Computing expectancies */ |
| hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); | |
| for(i=1; i<=nlstate;i++) | for(i=1; i<=nlstate;i++) |
| for(j=1; j<=nlstate;j++) | for(j=1; j<=nlstate;j++) |
| for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ | 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 ); | fprintf(ficresstdeij,"%3.0f",age ); |
| cptj=0; | 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(i=1; i<=nlstate;i++) |
| for(j=1; j<=nlstate;j++){ | for(j=1; j<=nlstate;j++){ |
| cptj++; | cptj= (j-1)*nlstate+i; |
| fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) ); | 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"); | printf("\n"); |
| fprintf(ficlog,"\n"); | fprintf(ficlog,"\n"); |
| free_vector(xm,1,npar); | |
| free_vector(xp,1,npar); | free_vector(xp,1,npar); |
| free_matrix(dnewm,1,nlstate*nlstate,1,npar); | free_matrix(dnewm,1,nlstate*nlstate,1,npar); |
| free_matrix(doldm,1,nlstate*nlstate,1,nlstate*nlstate); | free_matrix(doldm,1,nlstate*nlstate,1,nlstate*nlstate); |
| Line 2444 void evsij(char fileres[], double ***eij | Line 2709 void evsij(char fileres[], double ***eij |
| } | } |
| /************ Variance ******************/ | /************ 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 */ | /* Variance of health expectancies */ |
| /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ | /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ |
| Line 2495 void varevsij(char optionfilefiname[], d | Line 2760 void varevsij(char optionfilefiname[], d |
| fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev); | fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev); |
| } | } |
| printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%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); | 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,"# 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); | fprintf(ficresprobmorprev,"# Age cov=%-d",ij); |
| for(j=nlstate+1; j<=(nlstate+ndeath);j++){ | for(j=nlstate+1; j<=(nlstate+ndeath);j++){ |
| Line 2505 void varevsij(char optionfilefiname[], d | Line 2772 void varevsij(char optionfilefiname[], d |
| } | } |
| fprintf(ficresprobmorprev,"\n"); | fprintf(ficresprobmorprev,"\n"); |
| fprintf(ficgp,"\n# Routine varevsij"); | fprintf(ficgp,"\n# Routine varevsij"); |
| /* fprintf(fichtm, "#Local time at start: %s", strstart);*/ | |
| fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n"); | fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n"); |
| fprintf(fichtm,"\n<br>%s <br>\n",digitp); | fprintf(fichtm,"\n<br>%s <br>\n",digitp); |
| /* } */ | /* } */ |
| varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); | varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); |
| pstamp(ficresvij); | |
| 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"); | 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"); | fprintf(ficresvij,"# Age"); |
| for(i=1; i<=nlstate;i++) | for(i=1; i<=nlstate;i++) |
| for(j=1; j<=nlstate;j++) | 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"); | fprintf(ficresvij,"\n"); |
| xp=vector(1,npar); | xp=vector(1,npar); |
| Line 2742 void varevsij(char optionfilefiname[], d | Line 3014 void varevsij(char optionfilefiname[], d |
| } /* end varevsij */ | } /* end varevsij */ |
| /************ Variance of prevlim ******************/ | /************ 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 */ | /* Variance of prevalence limit */ |
| /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ | /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ |
| Line 2755 void varprevlim(char fileres[], double * | Line 3027 void varprevlim(char fileres[], double * |
| double **gradg, **trgradg; | double **gradg, **trgradg; |
| double age,agelim; | double age,agelim; |
| int theta; | int theta; |
| fprintf(ficresvpl,"# Standard deviation of stable prevalences \n"); | pstamp(ficresvpl); |
| fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n"); | |
| fprintf(ficresvpl,"# Age"); | fprintf(ficresvpl,"# Age"); |
| for(i=1; i<=nlstate;i++) | for(i=1; i<=nlstate;i++) |
| fprintf(ficresvpl," %1d-%1d",i,i); | fprintf(ficresvpl," %1d-%1d",i,i); |
| Line 2825 void varprevlim(char fileres[], double * | Line 3098 void varprevlim(char fileres[], double * |
| } | } |
| /************ Variance of one-step probabilities ******************/ | /************ 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 i, j=0, i1, k1, l1, t, tj; |
| int k2, l2, j1, z1; | int k2, l2, j1, z1; |
| Line 2870 void varprob(char optionfilefiname[], do | Line 3143 void varprob(char optionfilefiname[], do |
| fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); | 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); | 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); | 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,"#One-step probabilities and stand. devi in ()\n"); |
| fprintf(ficresprob,"# Age"); | fprintf(ficresprob,"# Age"); |
| pstamp(ficresprobcov); | |
| fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); | fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); |
| fprintf(ficresprobcov,"# Age"); | fprintf(ficresprobcov,"# Age"); |
| pstamp(ficresprobcor); | |
| fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); | fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); |
| fprintf(ficresprobcov,"# Age"); | fprintf(ficresprobcor,"# Age"); |
| for(i=1; i<=nlstate;i++) | for(i=1; i<=nlstate;i++) |
| Line 3127 To be simple, these graphs help to under | Line 3402 To be simple, these graphs help to under |
| } | } |
| free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); | 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(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); | free_vector(xp,1,npar); |
| fclose(ficresprob); | fclose(ficresprob); |
| fclose(ficresprobcov); | fclose(ficresprobcov); |
| Line 3145 void printinghtml(char fileres[], char t | Line 3422 void printinghtml(char fileres[], char t |
| double jprev2, double mprev2,double anprev2){ | double jprev2, double mprev2,double anprev2){ |
| int jj1, k1, i1, cpt; | int jj1, k1, i1, cpt; |
| fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n \ | fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \ |
| <li><a href='#secondorder'>Result files (second order (variance)</a>\n \ | |
| </ul>"); | |
| fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n \ | |
| - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ", | - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ", |
| jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p")); | jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p")); |
| fprintf(fichtm,"\ | fprintf(fichtm,"\ |
| - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ", | - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ", |
| stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij")); | stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij")); |
| fprintf(fichtm,"\ | fprintf(fichtm,"\ |
| - Stable prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", | - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", |
| subdirf2(fileres,"pl"),subdirf2(fileres,"pl")); | subdirf2(fileres,"pl"),subdirf2(fileres,"pl")); |
| fprintf(fichtm,"\ | fprintf(fichtm,"\ |
| - Life expectancies by age and initial health status (estepm=%2d months): \ | - (a) Life expectancies by health status at initial age, (b) health expectancies by health status at initial age: ei., eij (estepm=%2d months): \ |
| <a href=\"%s\">%s</a> <br>\n</li>", | <a href=\"%s\">%s</a> <br>\n</li>", |
| estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e")); | estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e")); |
| fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>"); | fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>"); |
| m=cptcoveff; | m=cptcoveff; |
| Line 3175 fprintf(fichtm," \n<ul><li><b>Graphs</b> | Line 3456 fprintf(fichtm," \n<ul><li><b>Graphs</b> |
| fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); | fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
| } | } |
| /* Pij */ | /* Pij */ |
| fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s%d1.png<br> \ | fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s%d1.png\">%s%d1.png</a><br> \ |
| <img src=\"%s%d1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); | <img src=\"%s%d1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); |
| /* Quasi-incidences */ | /* Quasi-incidences */ |
| fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\ | fprintf(fichtm,"<br>- 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<br> \ | before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: <a href=\"%s%d2.png\">%s%d2.png</a><br> \ |
| <img src=\"%s%d2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); | <img src=\"%s%d2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); |
| /* Stable prevalence in each health state */ | /* Period (stable) prevalence in each health state */ |
| for(cpt=1; cpt<nlstate;cpt++){ | for(cpt=1; cpt<nlstate;cpt++){ |
| fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br> \ | fprintf(fichtm,"<br>- Period (stable) prevalence in each health state : <a href=\"%s%d%d.png\">%s%d%d.png</a><br> \ |
| <img src=\"%s%d%d.png\">",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1); | <img src=\"%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++) { | for(cpt=1; cpt<=nlstate;cpt++) { |
| fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): %s%d%d.png <br> \ | fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \ |
| <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1); | <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1); |
| } | } |
| } /* end i1 */ | } /* end i1 */ |
| }/* End k1 */ | }/* End k1 */ |
| Line 3196 fprintf(fichtm," \n<ul><li><b>Graphs</b> | Line 3477 fprintf(fichtm," \n<ul><li><b>Graphs</b> |
| fprintf(fichtm,"\ | fprintf(fichtm,"\ |
| \n<br><li><h4> Result files (second order: variances)</h4>\n\ | \n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\ |
| - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n", rfileres,rfileres); | - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n", rfileres,rfileres); |
| fprintf(fichtm," - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", | fprintf(fichtm," - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", |
| Line 3209 fprintf(fichtm," \n<ul><li><b>Graphs</b> | Line 3490 fprintf(fichtm," \n<ul><li><b>Graphs</b> |
| - Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", | - Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", |
| subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor")); | subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor")); |
| fprintf(fichtm,"\ | fprintf(fichtm,"\ |
| - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"%s\">%s</a><br>\n", | - Variances and covariances of health expectancies by age and <b>initial health status</b> (cov(e<sup>ij</sup>,e<sup>kl</sup>)(estepm=%2d months): \ |
| <a href=\"%s\">%s</a> <br>\n</li>", | |
| estepm,subdirf2(fileres,"cve"),subdirf2(fileres,"cve")); | |
| fprintf(fichtm,"\ | |
| - (a) Health expectancies by health status at initial age (e<sup>ij</sup>) and standard errors (in parentheses) (b) life expectancies and standard errors (e<sup>i.</sup>=e<sup>i1</sup>+e<sup>i2</sup>+...)(estepm=%2d months): \ | |
| <a href=\"%s\">%s</a> <br>\n</li>", | |
| 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): <a href=\"%s\">%s</a><br>\n", | |
| estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v")); | estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v")); |
| fprintf(fichtm,"\ | fprintf(fichtm,"\ |
| - Health expectancies with their variances (no covariance): <a href=\"%s\">%s</a> <br>\n", | - Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors: <a href=\"%s\">%s</a> <br>\n", |
| subdirf2(fileres,"t"),subdirf2(fileres,"t")); | subdirf2(fileres,"t"),subdirf2(fileres,"t")); |
| fprintf(fichtm,"\ | fprintf(fichtm,"\ |
| - Standard deviation of stable prevalences: <a href=\"%s\">%s</a> <br>\n",\ | - Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\ |
| subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl")); | subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl")); |
| /* if(popforecast==1) fprintf(fichtm,"\n */ | /* if(popforecast==1) fprintf(fichtm,"\n */ |
| Line 3287 plot [%.f:%.f] \"%s\" every :::%d::%d u | Line 3576 plot [%.f:%.f] \"%s\" every :::%d::%d u |
| if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); | if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); |
| else 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 ++) { | for (i=1; i<= nlstate ; i ++) { |
| if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); | if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); |
| else fprintf(ficgp," \%%*lf (\%%*lf)"); | else fprintf(ficgp," \%%*lf (\%%*lf)"); |
| Line 3335 plot [%.f:%.f] \"%s\" every :::%d::%d u | Line 3624 plot [%.f:%.f] \"%s\" every :::%d::%d u |
| for (k1=1; k1<= m ; k1 ++) { | for (k1=1; k1<= m ; k1 ++) { |
| for (cpt=1; cpt<= nlstate ; cpt ++) { | 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,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1); |
| fprintf(ficgp,"set ter png small\n\ | fprintf(ficgp,"set ter png small\n\ |
| set size 0.65,0.65\n\ | set size 0.65,0.65\n\ |
| Line 3349 plot [%.f:%.f] \"%s\" every :::%d::%d u | Line 3639 plot [%.f:%.f] \"%s\" every :::%d::%d u |
| */ | */ |
| for (i=1; i< nlstate ; i ++) { | 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); | |
| } | } |
| } | } |
| Line 3910 double gompertz(double x[]) | Line 4202 double gompertz(double x[]) |
| { | { |
| double A,B,L=0.0,sump=0.,num=0.; | double A,B,L=0.0,sump=0.,num=0.; |
| int i,n=0; /* n is the size of the sample */ | int i,n=0; /* n is the size of the sample */ |
| for (i=0;i<=imx-1 ; i++) { | for (i=0;i<=imx-1 ; i++) { |
| sump=sump+weight[i]; | sump=sump+weight[i]; |
| sump=sump+1; | /* sump=sump+1;*/ |
| num=num+1; | num=num+1; |
| } | } |
| /* for (i=1; i<=imx; i++) | /* 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]);*/ | if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/ |
| for (i=0;i<=imx-1 ; i++) | for (i=1;i<=imx ; i++) |
| { | { |
| if (cens[i]==1 & wav[i]>1) | if (cens[i] == 1 && wav[i]>1) |
| A=-x[1]/(x[2])* | A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); |
| (exp(x[2]/YEARM*(agecens[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12))); | |
| if (cens[i]==0 & wav[i]>1) | if (cens[i] == 0 && wav[i]>1) |
| A=-x[1]/(x[2])* | A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))) |
| (exp(x[2]/YEARM*(agedc[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12))) | +log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM); |
| +log(x[1]/YEARM)+x[2]/YEARM*(agedc[i]*12-agegomp*12)+log(YEARM); | |
| if (wav[i]>1 & agecens[i]>15) { | /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */ |
| if (wav[i] > 1 ) { /* ??? */ | |
| L=L+A*weight[i]; | 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("\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]);*/ |
| } | } |
| Line 3945 double gompertz(double x[]) | Line 4237 double gompertz(double x[]) |
| /******************* Printing html file ***********/ | /******************* Printing html file ***********/ |
| void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \ | void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \ |
| int lastpass, int stepm, int weightopt, char model[],\ | int lastpass, int stepm, int weightopt, char model[],\ |
| int imx, double p[],double **matcov){ | int imx, double p[],double **matcov,double agemortsup){ |
| int i; | int i,k; |
| fprintf(fichtm,"<ul><li><h4>Result files </h4>\n Force of mortality. Parameters of the Gompertz fit (with confidence interval in brackets):<br>"); | fprintf(fichtm,"<ul><li><h4>Result files </h4>\n Force of mortality. Parameters of the Gompertz fit (with confidence interval in brackets):<br>"); |
| fprintf(fichtm," mu(age) =%lf*exp(%lf*(age-%d)) per year<br><br>",p[1],p[2],agegomp); | fprintf(fichtm," mu(age) =%lf*exp(%lf*(age-%d)) per year<br><br>",p[1],p[2],agegomp); |
| Line 3954 void printinghtmlmort(char fileres[], ch | Line 4246 void printinghtmlmort(char fileres[], ch |
| fprintf(fichtm," p[%d] = %lf [%f ; %f]<br>\n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); | fprintf(fichtm," p[%d] = %lf [%f ; %f]<br>\n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); |
| fprintf(fichtm,"<br><br><img src=\"graphmort.png\">"); | fprintf(fichtm,"<br><br><img src=\"graphmort.png\">"); |
| fprintf(fichtm,"</ul>"); | fprintf(fichtm,"</ul>"); |
| fprintf(fichtm,"<ul><li><h4>Life table</h4>\n <br>"); | |
| fprintf(fichtm,"\nAge l<inf>x</inf> q<inf>x</inf> d(x,x+1) L<inf>x</inf> T<inf>x</inf> e<infx</inf><br>"); | |
| for (k=agegomp;k<(agemortsup-2);k++) | |
| fprintf(fichtm,"%d %.0lf %lf %.0lf %.0lf %.0lf %lf<br>\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); | fflush(fichtm); |
| } | } |
| Line 3983 void printinggnuplotmort(char fileres[], | Line 4284 void printinggnuplotmort(char fileres[], |
| /***********************************************/ | /***********************************************/ |
| /**************** Main Program *****************/ | /**************** Main Program *****************/ |
| /***********************************************/ | /***********************************************/ |
| Line 3991 int main(int argc, char *argv[]) | Line 4293 int main(int argc, char *argv[]) |
| { | { |
| int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav); | 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 i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod; |
| int linei, month, year,iout; | |
| int jj, ll, li, lj, lk, imk; | int jj, ll, li, lj, lk, imk; |
| int numlinepar=0; /* Current linenumber of parameter file */ | int numlinepar=0; /* Current linenumber of parameter file */ |
| int itimes; | int itimes; |
| int NDIM=2; | int NDIM=2; |
| char ca[32], cb[32], cc[32]; | char ca[32], cb[32], cc[32]; |
| char dummy[]=" "; | |
| /* FILE *fichtm; *//* Html File */ | /* FILE *fichtm; *//* Html File */ |
| /* FILE *ficgp;*/ /*Gnuplot File */ | /* FILE *ficgp;*/ /*Gnuplot File */ |
| struct stat info; | |
| double agedeb, agefin,hf; | double agedeb, agefin,hf; |
| double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; | double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; |
| Line 4012 int main(int argc, char *argv[]) | Line 4317 int main(int argc, char *argv[]) |
| char line[MAXLINE], linepar[MAXLINE]; | char line[MAXLINE], linepar[MAXLINE]; |
| char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; | char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; |
| char pathr[MAXLINE], pathimach[MAXLINE]; | char pathr[MAXLINE], pathimach[MAXLINE]; |
| char **bp, *tok, *val; /* pathtot */ | |
| int firstobs=1, lastobs=10; | int firstobs=1, lastobs=10; |
| int sdeb, sfin; /* Status at beginning and end */ | int sdeb, sfin; /* Status at beginning and end */ |
| int c, h , cpt,l; | int c, h , cpt,l; |
| Line 4021 int main(int argc, char *argv[]) | Line 4327 int main(int argc, char *argv[]) |
| int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ | int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ |
| int mobilav=0,popforecast=0; | int mobilav=0,popforecast=0; |
| int hstepm, nhstepm; | int hstepm, nhstepm; |
| int agemortsup; | |
| float sumlpop=0.; | |
| double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000; | 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; | double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000; |
| Line 4045 int main(int argc, char *argv[]) | Line 4353 int main(int argc, char *argv[]) |
| char z[1]="c", occ; | char z[1]="c", occ; |
| char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80]; | 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; | char *stratrunc; |
| int lstra; | int lstra; |
| Line 4087 int main(int argc, char *argv[]) | Line 4395 int main(int argc, char *argv[]) |
| printf("\n%s\n%s",version,fullversion); | printf("\n%s\n%s",version,fullversion); |
| if(argc <=1){ | if(argc <=1){ |
| printf("\nEnter the parameter file name: "); | 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; /* Un peu sale */ | |
| } | |
| } | } |
| else{ | else{ |
| strcpy(pathtot,argv[1]); | strcpy(pathtot,argv[1]); |
| Line 4097 int main(int argc, char *argv[]) | Line 4415 int main(int argc, char *argv[]) |
| printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/ | printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/ |
| /* cutv(path,optionfile,pathtot,'\\');*/ | /* 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); | 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]); */ | /* strcpy(pathimach,argv[0]); */ |
| /* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */ | |
| split(pathtot,path,optionfile,optionfilext,optionfilefiname); | split(pathtot,path,optionfile,optionfilext,optionfilefiname); |
| printf("pathimach=%s, pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname); | printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); |
| chdir(path); | chdir(path); |
| strcpy(command,"mkdir "); | strcpy(command,"mkdir "); |
| strcat(command,optionfilefiname); | strcat(command,optionfilefiname); |
| Line 4267 int main(int argc, char *argv[]) | Line 4589 int main(int argc, char *argv[]) |
| } | } |
| fflush(ficlog); | fflush(ficlog); |
| p=param[1][1]; | p=param[1][1]; |
| /* Reads comments: lines beginning with '#' */ | /* Reads comments: lines beginning with '#' */ |
| Line 4391 int main(int argc, char *argv[]) | Line 4712 int main(int argc, char *argv[]) |
| ncodemax=ivector(1,8); | ncodemax=ivector(1,8); |
| i=1; | i=1; |
| while (fgets(line, MAXLINE, fic) != NULL) { | linei=0; |
| if ((i >= firstobs) && (i <=lastobs)) { | while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) { |
| linei=linei+1; | |
| for (j=maxwav;j>=1;j--){ | for(j=strlen(line); j>=0;j--){ /* Untabifies line */ |
| cutv(stra, strb,line,' '); s[j][i]=atoi(strb); | if(line[j] == '\t') |
| strcpy(line,stra); | line[j] = ' '; |
| 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); | 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); | |
| } | } |
| s[j][i]=lval; | |
| cutv(stra, strb,line,'/'); andc[i]=(double)(atoi(strb)); strcpy(line,stra); | |
| cutv(stra, strb,line,' '); moisdc[i]=(double)(atoi(strb)); strcpy(line,stra); | strcpy(line,stra); |
| cutv(stra, strb,line,' '); | |
| cutv(stra, strb,line,'/'); annais[i]=(double)(atoi(strb)); strcpy(line,stra); | if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ |
| cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra); | } |
| else if(iout=sscanf(strb,"%s.") != 0){ | |
| cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra); | month=99; |
| for (j=ncovcol;j>=1;j--){ | year=9999; |
| cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra); | }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); |
| lstra=strlen(stra); | exit(1); |
| if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */ | } |
| stratrunc = &(stra[lstra-9]); | anint[j][i]= (double) year; |
| num[i]=atol(stratrunc); | 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); | |
| } | } |
| else | covar[j][i]=(double)(lval); |
| num[i]=atol(stra); | strcpy(line,stra); |
| } | |
| /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){ | lstra=strlen(stra); |
| 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;}*/ | |
| if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */ | |
| i=i+1; | 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); | /* printf("ii=%d", ij); |
| scanf("%d",i);*/ | scanf("%d",i);*/ |
| imx=i-1; /* Number of individuals */ | imx=i-1; /* Number of individuals */ |
| Line 4438 int main(int argc, char *argv[]) | Line 4836 int main(int argc, char *argv[]) |
| if (s[4][i]==9) s[4][i]=-1; | 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]));}*/ | 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; | /*if ((s[3][i]==3) || (s[4][i]==3)) weight[i]=0.08; |
| else weight[i]=1;*/ | 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 */ | Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */ |
| Tprod=ivector(1,15); | Tprod=ivector(1,15); |
| Tvaraff=ivector(1,15); | Tvaraff=ivector(1,15); |
| Line 4522 int main(int argc, char *argv[]) | Line 4920 int main(int argc, char *argv[]) |
| printf("cptcovprod=%d ", cptcovprod); | printf("cptcovprod=%d ", cptcovprod); |
| fprintf(ficlog,"cptcovprod=%d ", cptcovprod); | fprintf(ficlog,"cptcovprod=%d ", cptcovprod); |
| scanf("%d ",i); | scanf("%d ",i);*/ |
| fclose(fic);*/ | |
| /* if(mle==1){*/ | /* if(mle==1){*/ |
| if (weightopt != 1) { /* Maximisation without weights*/ | if (weightopt != 1) { /* Maximisation without weights*/ |
| Line 4556 int main(int argc, char *argv[]) | Line 4953 int main(int argc, char *argv[]) |
| for (i=1; i<=imx; i++) { | for (i=1; i<=imx; i++) { |
| agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); | agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); |
| for(m=firstpass; (m<= lastpass); m++){ | 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 (s[m][i] >= nlstate+1) { |
| if(agedc[i]>0) | if(agedc[i]>0) |
| if((int)moisdc[i]!=99 && (int)andc[i]!=9999) | if((int)moisdc[i]!=99 && (int)andc[i]!=9999) |
| Line 4572 int main(int argc, char *argv[]) | Line 4969 int main(int argc, char *argv[]) |
| } | } |
| } | } |
| else if(s[m][i] !=9){ /* Standard case, age in fractional | else if(s[m][i] !=9){ /* Standard case, age in fractional |
| years but with the precision of a | years but with the precision of a month */ |
| month */ | |
| agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]); | 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) | if((int)mint[m][i]==99 || (int)anint[m][i]==9999) |
| agev[m][i]=1; | agev[m][i]=1; |
| Line 4715 Title=%s <br>Datafile=%s Firstpass=%d La | Line 5111 Title=%s <br>Datafile=%s Firstpass=%d La |
| \n\ | \n\ |
| <hr size=\"2\" color=\"#EC5E5E\">\ | <hr size=\"2\" color=\"#EC5E5E\">\ |
| <ul><li><h4>Parameter files</h4>\n\ | <ul><li><h4>Parameter files</h4>\n\ |
| - Parameter file: <a href=\"%s.%s\">%s.%s</a><br>\n\ | |
| - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n\ | - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n\ |
| - Log file of the run: <a href=\"%s\">%s</a><br>\n\ | - Log file of the run: <a href=\"%s\">%s</a><br>\n\ |
| - Gnuplot file name: <a href=\"%s\">%s</a><br>\n\ | - Gnuplot file name: <a href=\"%s\">%s</a><br>\n\ |
| - Date and time at start: %s</ul>\n",\ | - Date and time at start: %s</ul>\n",\ |
| fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\ | fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\ |
| optionfilefiname,optionfilext,optionfilefiname,optionfilext,\ | |
| fileres,fileres,\ | fileres,fileres,\ |
| filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart); | filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart); |
| fflush(fichtm); | fflush(fichtm); |
| Line 4730 Title=%s <br>Datafile=%s Firstpass=%d La | Line 5128 Title=%s <br>Datafile=%s Firstpass=%d La |
| /* Calculates basic frequencies. Computes observed prevalence at single age | /* Calculates basic frequencies. Computes observed prevalence at single age |
| and prints on file fileres'p'. */ | 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,"\n"); |
| fprintf(fichtm,"<br>Total number of observations=%d <br>\n\ | fprintf(fichtm,"<br>Total number of observations=%d <br>\n\ |
| Line 4749 Interval (in months) between two waves: | Line 5147 Interval (in months) between two waves: |
| p=param[1][1]; /* *(*(*(param +1)+1)+0) */ | p=param[1][1]; /* *(*(*(param +1)+1)+0) */ |
| globpr=0; /* To get the number ipmx of contributions and the sum of weights*/ | globpr=0; /* To get the number ipmx of contributions and the sum of weights*/ |
| if (mle==-3){ | if (mle==-3){ |
| ximort=matrix(1,NDIM,1,NDIM); | ximort=matrix(1,NDIM,1,NDIM); |
| cens=ivector(1,n); | cens=ivector(1,n); |
| Line 4758 Interval (in months) between two waves: | Line 5157 Interval (in months) between two waves: |
| for (i=1; i<=imx; i++){ | for (i=1; i<=imx; i++){ |
| dcwave[i]=-1; | dcwave[i]=-1; |
| for (j=1; j<=lastpass; j++) | for (m=firstpass; m<=lastpass; m++) |
| if (s[j][i]>nlstate) { | if (s[m][i]>nlstate) { |
| dcwave[i]=j; | dcwave[i]=m; |
| /* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/ | /* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/ |
| break; | break; |
| } | } |
| Line 4769 Interval (in months) between two waves: | Line 5168 Interval (in months) between two waves: |
| for (i=1; i<=imx; i++) { | for (i=1; i<=imx; i++) { |
| if (wav[i]>0){ | if (wav[i]>0){ |
| ageexmed[i]=agev[mw[1][i]][i]; | ageexmed[i]=agev[mw[1][i]][i]; |
| j=wav[i];agecens[i]=1.; | j=wav[i]; |
| if (ageexmed[i]>1 & wav[i]>0) agecens[i]=agev[mw[j][i]][i]; | agecens[i]=1.; |
| cens[i]=1; | |
| if (ageexmed[i]> 1 && wav[i] > 0){ | |
| if (ageexmed[i]<1) cens[i]=-1; | agecens[i]=agev[mw[j][i]][i]; |
| if (agedc[i]< AGESUP & agedc[i]>1 & dcwave[i]>firstpass & dcwave[i]<=lastpass) cens[i]=0 ; | 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; | else cens[i]=-1; |
| } | } |
| Line 4783 Interval (in months) between two waves: | Line 5186 Interval (in months) between two waves: |
| for (j=1;j<=NDIM;j++) | for (j=1;j<=NDIM;j++) |
| ximort[i][j]=(i == j ? 1.0 : 0.0); | ximort[i][j]=(i == j ? 1.0 : 0.0); |
| } | } |
| p[1]=0.1; p[2]=0.1; | p[1]=0.0268; p[NDIM]=0.083; |
| /*printf("%lf %lf", p[1], p[2]);*/ | /*printf("%lf %lf", p[1], p[2]);*/ |
| printf("Powell\n"); fprintf(ficlog,"Powell\n"); | printf("Powell\n"); fprintf(ficlog,"Powell\n"); |
| strcpy(filerespow,"pow-mort"); | strcpy(filerespow,"pow-mort"); |
| strcat(filerespow,fileres); | strcat(filerespow,fileres); |
| if((ficrespow=fopen(filerespow,"w"))==NULL) { | if((ficrespow=fopen(filerespow,"w"))==NULL) { |
| printf("Problem with resultfile: %s\n", filerespow); | printf("Problem with resultfile: %s\n", filerespow); |
| fprintf(ficlog,"Problem with resultfile: %s\n", filerespow); | fprintf(ficlog,"Problem with resultfile: %s\n", filerespow); |
| } | } |
| fprintf(ficrespow,"# Powell\n# iter -2*LL"); | fprintf(ficrespow,"# Powell\n# iter -2*LL"); |
| /* for (i=1;i<=nlstate;i++) | /* for (i=1;i<=nlstate;i++) |
| for(j=1;j<=nlstate+ndeath;j++) | for(j=1;j<=nlstate+ndeath;j++) |
| if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); | if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); |
| */ | */ |
| fprintf(ficrespow,"\n"); | fprintf(ficrespow,"\n"); |
| powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); | powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); |
| fclose(ficrespow); | fclose(ficrespow); |
| hesscov(matcov, p, NDIM,delti, 1e-4, gompertz); | hesscov(matcov, p, NDIM, delti, 1e-4, gompertz); |
| for(i=1; i <=NDIM; i++) | for(i=1; i <=NDIM; i++) |
| for(j=i+1;j<=NDIM;j++) | for(j=i+1;j<=NDIM;j++) |
| Line 4822 Interval (in months) between two waves: | Line 5225 Interval (in months) between two waves: |
| printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp); | 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++) | 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])); | 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<agemortsup;k++) | |
| lsurv[k+1]=lsurv[k]-lsurv[k]*(p[1]*exp(p[2]*(k-agegomp))); | |
| for (k=agegomp;k<agemortsup;k++){ | |
| lpop[k]=(lsurv[k]+lsurv[k+1])/2.; | |
| sumlpop=sumlpop+lpop[k]; | |
| } | |
| tpop[agegomp]=sumlpop; | |
| for (k=agegomp;k<(agemortsup-3);k++){ | |
| /* tpop[k+1]=2;*/ | |
| tpop[k+1]=tpop[k]-lpop[k]; | |
| } | |
| printf("\nAge lx qx dx Lx Tx e(x)\n"); | |
| for (k=agegomp;k<(agemortsup-2);k++) | |
| printf("%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]); | |
| replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */ | replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */ |
| printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); | printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); |
| printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \ | printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \ |
| stepm, weightopt,\ | stepm, weightopt,\ |
| model,imx,p,matcov); | model,imx,p,matcov,agemortsup); |
| free_vector(lsurv,1,AGESUP); | |
| free_vector(lpop,1,AGESUP); | |
| free_vector(tpop,1,AGESUP); | |
| } /* Endof if mle==-3 */ | } /* Endof if mle==-3 */ |
| else{ /* For mle >=1 */ | else{ /* For mle >=1 */ |
| likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ | likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ |
| Line 5077 Interval (in months) between two waves: | Line 5515 Interval (in months) between two waves: |
| fclose(ficres); | fclose(ficres); |
| /*--------------- Prevalence limit (stable prevalence) --------------*/ | /*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
| strcpy(filerespl,"pl"); | strcpy(filerespl,"pl"); |
| strcat(filerespl,fileres); | strcat(filerespl,fileres); |
| if((ficrespl=fopen(filerespl,"w"))==NULL) { | if((ficrespl=fopen(filerespl,"w"))==NULL) { |
| printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end; | printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end; |
| fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end; | fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end; |
| } | } |
| printf("Computing stable prevalence: result on file '%s' \n", filerespl); | printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl); |
| fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl); | fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl); |
| fprintf(ficrespl,"#Stable prevalence \n"); | pstamp(ficrespl); |
| fprintf(ficrespl,"# Period (stable) prevalence \n"); | |
| fprintf(ficrespl,"#Age "); | fprintf(ficrespl,"#Age "); |
| for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); | for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); |
| fprintf(ficrespl,"\n"); | fprintf(ficrespl,"\n"); |
| Line 5147 Interval (in months) between two waves: | Line 5586 Interval (in months) between two waves: |
| hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ | hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ |
| /* hstepm=1; aff par mois*/ | /* 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 "); | 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(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ | for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
| Line 5184 Interval (in months) between two waves: | Line 5623 Interval (in months) between two waves: |
| } | } |
| } | } |
| 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); |
| Line 5217 Interval (in months) between two waves: | Line 5656 Interval (in months) between two waves: |
| printf("Problem with total LE resultfile: %s\n", filerest);goto end; | printf("Problem with total LE resultfile: %s\n", filerest);goto end; |
| fprintf(ficlog,"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); | printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); |
| fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest); | fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); |
| strcpy(filerese,"e"); | strcpy(filerese,"e"); |
| Line 5230 Interval (in months) between two waves: | Line 5669 Interval (in months) between two waves: |
| printf("Computing Health Expectancies: result on file '%s' \n", filerese); | printf("Computing Health Expectancies: result on file '%s' \n", filerese); |
| fprintf(ficlog,"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"); | strcpy(fileresv,"v"); |
| strcat(fileresv,fileres); | strcat(fileresv,fileres); |
| if((ficresvij=fopen(fileresv,"w"))==NULL) { | if((ficresvij=fopen(fileresv,"w"))==NULL) { |
| Line 5262 Interval (in months) between two waves: | Line 5719 Interval (in months) between two waves: |
| fprintf(ficrest,"******\n"); | fprintf(ficrest,"******\n"); |
| fprintf(ficreseij,"\n#****** "); | fprintf(ficreseij,"\n#****** "); |
| for(j=1;j<=cptcoveff;j++) | 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(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(ficreseij,"******\n"); |
| fprintf(ficresstdeij,"******\n"); | |
| fprintf(ficrescveij,"******\n"); | |
| fprintf(ficresvij,"\n#****** "); | fprintf(ficresvij,"\n#****** "); |
| for(j=1;j<=cptcoveff;j++) | for(j=1;j<=cptcoveff;j++) |
| Line 5273 Interval (in months) between two waves: | Line 5737 Interval (in months) between two waves: |
| eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); | eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); |
| oldm=oldms;savm=savms; | oldm=oldms;savm=savms; |
| evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov); | 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); | |
| vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); | vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); |
| oldm=oldms;savm=savms; | 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); | 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){ | 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); | varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav, strstart); |
| } | } |
| pstamp(ficrest); | |
| fprintf(ficrest,"#Total LEs with variances: e.. (std) "); | fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state ( e.. (std) "); |
| for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); | for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); |
| fprintf(ficrest,"\n"); | fprintf(ficrest,"\n"); |
| Line 5331 Interval (in months) between two waves: | Line 5796 Interval (in months) between two waves: |
| free_ivector(cod,1,n); | free_ivector(cod,1,n); |
| free_ivector(tab,1,NCOVMAX); | free_ivector(tab,1,NCOVMAX); |
| fclose(ficreseij); | fclose(ficreseij); |
| fclose(ficresstdeij); | |
| fclose(ficrescveij); | |
| fclose(ficresvij); | fclose(ficresvij); |
| fclose(ficrest); | fclose(ficrest); |
| fclose(ficpar); | fclose(ficpar); |
| /*------- Variance of stable prevalence------*/ | /*------- Variance of period (stable) prevalence------*/ |
| strcpy(fileresvpl,"vpl"); | strcpy(fileresvpl,"vpl"); |
| strcat(fileresvpl,fileres); | strcat(fileresvpl,fileres); |
| if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { | if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { |
| printf("Problem with variance of stable prevalence resultfile: %s\n", fileresvpl); | printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl); |
| exit(0); | exit(0); |
| } | } |
| printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl); | printf("Computing Variance-covariance of period (stable) prevalence: file '%s' \n", fileresvpl); |
| for(cptcov=1,k=0;cptcov<=i1;cptcov++){ | for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ | for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
| Line 5355 Interval (in months) between two waves: | Line 5822 Interval (in months) between two waves: |
| varpl=matrix(1,nlstate,(int) bage, (int) fage); | varpl=matrix(1,nlstate,(int) bage, (int) fage); |
| oldm=oldms;savm=savms; | oldm=oldms;savm=savms; |
| varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k); | 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); | free_matrix(varpl,1,nlstate,(int) bage, (int)fage); |
| } | } |
| } | } |
| Line 5367 Interval (in months) between two waves: | Line 5834 Interval (in months) between two waves: |
| free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); | free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| } /* mle==-3 arrives here for freeing */ | } /* mle==-3 arrives here for freeing */ |
| free_matrix(prlim,1,nlstate,1,nlstate); | |
| free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); | free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(oldms, 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(newms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); | free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(covar,0,NCOVMAX,1,n); | free_matrix(covar,0,NCOVMAX,1,n); |
| free_matrix(matcov,1,npar,1,npar); | free_matrix(matcov,1,npar,1,npar); |
| /*free_vector(delti,1,npar);*/ | /*free_vector(delti,1,npar);*/ |
| Line 5386 Interval (in months) between two waves: | Line 5853 Interval (in months) between two waves: |
| free_ivector(Tage,1,15); | free_ivector(Tage,1,15); |
| free_ivector(Tcode,1,100); | free_ivector(Tcode,1,100); |
| free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX); | |
| free_imatrix(codtab,1,100,1,10); | |
| fflush(fichtm); | fflush(fichtm); |
| fflush(ficgp); | fflush(ficgp); |
| Line 5404 Interval (in months) between two waves: | Line 5872 Interval (in months) between two waves: |
| tm = *localtime(&end_time.tv_sec); | tm = *localtime(&end_time.tv_sec); |
| tmg = *gmtime(&end_time.tv_sec); | tmg = *gmtime(&end_time.tv_sec); |
| strcpy(strtend,asctime(&tm)); | strcpy(strtend,asctime(&tm)); |
| printf("Local time at start %s\nLocaltime at end %s",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); | 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 used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout)); |
| Line 5421 Interval (in months) between two waves: | Line 5889 Interval (in months) between two waves: |
| /*------ End -----------*/ | /*------ End -----------*/ |
| chdir(path); | chdir(path); |
| strcpy(plotcmd,"\""); | /*strcat(plotcmd,CHARSEPARATOR);*/ |
| strcat(plotcmd,pathimach); | sprintf(plotcmd,"gnuplot"); |
| strcat(plotcmd,GNUPLOTPROGRAM); | #ifndef UNIX |
| strcat(plotcmd,"\""); | sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach); |
| strcat(plotcmd," "); | #endif |
| strcat(plotcmd,optionfilegnuplot); | if(!stat(plotcmd,&info)){ |
| printf("Starting graphs with: %s",plotcmd);fflush(stdout); | 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){ | if((outcmd=system(plotcmd)) != 0){ |
| printf(" Problem with gnuplot\n"); | printf("\n Problem with gnuplot\n"); |
| } | } |
| printf(" Wait..."); | printf(" Wait..."); |
| while (z[0] != 'q') { | while (z[0] != 'q') { |