|
|
| version 1.91, 2003/06/25 15:30:29 | version 1.105, 2006/01/05 20:23:19 |
|---|---|
| Line 1 | Line 1 |
| /* $Id$ | /* $Id$ |
| $State$ | $State$ |
| $Log$ | $Log$ |
| Revision 1.105 2006/01/05 20:23:19 lievre | |
| *** empty log message *** | |
| Revision 1.104 2005/09/30 16:11:43 lievre | |
| (Module): sump fixed, loop imx fixed, and simplifications. | |
| (Module): If the status is missing at the last wave but we know | |
| that the person is alive, then we can code his/her status as -2 | |
| (instead of missing=-1 in earlier versions) and his/her | |
| contributions to the likelihood is 1 - Prob of dying from last | |
| health status (= 1-p13= p11+p12 in the easiest case of somebody in | |
| the healthy state at last known wave). Version is 0.98 | |
| Revision 1.103 2005/09/30 15:54:49 lievre | |
| (Module): sump fixed, loop imx fixed, and simplifications. | |
| Revision 1.102 2004/09/15 17:31:30 brouard | |
| Add the possibility to read data file including tab characters. | |
| Revision 1.101 2004/09/15 10:38:38 brouard | |
| Fix on curr_time | |
| Revision 1.100 2004/07/12 18:29:06 brouard | |
| Add version for Mac OS X. Just define UNIX in Makefile | |
| Revision 1.99 2004/06/05 08:57:40 brouard | |
| *** empty log message *** | |
| Revision 1.98 2004/05/16 15:05:56 brouard | |
| New version 0.97 . First attempt to estimate force of mortality | |
| directly from the data i.e. without the need of knowing the health | |
| state at each age, but using a Gompertz model: log u =a + b*age . | |
| This is the basic analysis of mortality and should be done before any | |
| other analysis, in order to test if the mortality estimated from the | |
| cross-longitudinal survey is different from the mortality estimated | |
| from other sources like vital statistic data. | |
| The same imach parameter file can be used but the option for mle should be -3. | |
| Agnès, who wrote this part of the code, tried to keep most of the | |
| former routines in order to include the new code within the former code. | |
| The output is very simple: only an estimate of the intercept and of | |
| the slope with 95% confident intervals. | |
| Current limitations: | |
| A) Even if you enter covariates, i.e. with the | |
| model= V1+V2 equation for example, the programm does only estimate a unique global model without covariates. | |
| B) There is no computation of Life Expectancy nor Life Table. | |
| Revision 1.97 2004/02/20 13:25:42 lievre | |
| Version 0.96d. Population forecasting command line is (temporarily) | |
| suppressed. | |
| Revision 1.96 2003/07/15 15:38:55 brouard | |
| * imach.c (Repository): Errors in subdirf, 2, 3 while printing tmpout is | |
| rewritten within the same printf. Workaround: many printfs. | |
| Revision 1.95 2003/07/08 07:54:34 brouard | |
| * imach.c (Repository): | |
| (Repository): Using imachwizard code to output a more meaningful covariance | |
| matrix (cov(a12,c31) instead of numbers. | |
| Revision 1.94 2003/06/27 13:00:02 brouard | |
| Just cleaning | |
| Revision 1.93 2003/06/25 16:33:55 brouard | |
| (Module): On windows (cygwin) function asctime_r doesn't | |
| exist so I changed back to asctime which exists. | |
| (Module): Version 0.96b | |
| Revision 1.92 2003/06/25 16:30:45 brouard | |
| (Module): On windows (cygwin) function asctime_r doesn't | |
| exist so I changed back to asctime which exists. | |
| Revision 1.91 2003/06/25 15:30:29 brouard | Revision 1.91 2003/06/25 15:30:29 brouard |
| * imach.c (Repository): Duplicated warning errors corrected. | * imach.c (Repository): Duplicated warning errors corrected. |
| (Repository): Elapsed time after each iteration is now output. It | (Repository): Elapsed time after each iteration is now output. It |
| Line 149 | Line 223 |
| #include <stdlib.h> | #include <stdlib.h> |
| #include <unistd.h> | #include <unistd.h> |
| #include <sys/time.h> | /* #include <sys/time.h> */ |
| #include <time.h> | #include <time.h> |
| #include "timeval.h" | #include "timeval.h" |
| /* #include <libintl.h> */ | |
| /* #define _(String) gettext (String) */ | |
| #define MAXLINE 256 | #define MAXLINE 256 |
| #define GNUPLOTPROGRAM "gnuplot" | #define GNUPLOTPROGRAM "gnuplot" |
| /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ | /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ |
| Line 173 | Line 250 |
| #define YEARM 12. /* Number of months per year */ | #define YEARM 12. /* Number of months per year */ |
| #define AGESUP 130 | #define AGESUP 130 |
| #define AGEBASE 40 | #define AGEBASE 40 |
| #ifdef unix | #define AGEGOMP 10. /* Minimal age for Gompertz adjustment */ |
| #ifdef UNIX | |
| #define DIRSEPARATOR '/' | #define DIRSEPARATOR '/' |
| #define ODIRSEPARATOR '\\' | #define ODIRSEPARATOR '\\' |
| #else | #else |
| Line 184 | Line 262 |
| /* $Id$ */ | /* $Id$ */ |
| /* $State$ */ | /* $State$ */ |
| char version[]="Imach version 0.96a, June 2003, INED-EUROREVES "; | char version[]="Imach version 0.98, September 2005, INED-EUROREVES "; |
| char fullversion[]="$Revision$ $Date$"; | char fullversion[]="$Revision$ $Date$"; |
| 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; |
| Line 214 int globpr; /* Global variable for print | Line 292 int globpr; /* Global variable for print |
| double fretone; /* Only one call to likelihood */ | double fretone; /* Only one call to likelihood */ |
| long ipmx; /* Number of contributions */ | long ipmx; /* Number of contributions */ |
| double sw; /* Sum of weights */ | double sw; /* Sum of weights */ |
| char filerespow[FILENAMELENGTH]; | |
| char fileresilk[FILENAMELENGTH]; /* File of individual contributions to the likelihood */ | char fileresilk[FILENAMELENGTH]; /* File of individual contributions to the likelihood */ |
| FILE *ficresilk; | FILE *ficresilk; |
| FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor; | FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor; |
| Line 228 char fileresvpl[FILENAMELENGTH]; | Line 307 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 optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH]; |
| char tmpout[FILENAMELENGTH]; | char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH]; |
| char command[FILENAMELENGTH]; | char command[FILENAMELENGTH]; |
| int outcmd=0; | int outcmd=0; |
| char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH]; | char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH]; |
| char lfileres[FILENAMELENGTH]; | |
| char filelog[FILENAMELENGTH]; /* Log file */ | char filelog[FILENAMELENGTH]; /* Log file */ |
| char filerest[FILENAMELENGTH]; | char filerest[FILENAMELENGTH]; |
| char fileregp[FILENAMELENGTH]; | char fileregp[FILENAMELENGTH]; |
| Line 276 static double maxarg1,maxarg2; | Line 355 static double maxarg1,maxarg2; |
| static double sqrarg; | static double sqrarg; |
| #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 :sqrarg*sqrarg) | #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 :sqrarg*sqrarg) |
| #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;} | #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;} |
| int agegomp= AGEGOMP; | |
| int imx; | int imx; |
| int stepm; | int stepm=1; |
| /* Stepm, step in month: minimum step interpolation*/ | /* Stepm, step in month: minimum step interpolation*/ |
| int estepm; | int estepm; |
| Line 286 int estepm; | Line 366 int estepm; |
| int m,nb; | int m,nb; |
| long *num; | long *num; |
| int firstpass=0, lastpass=4,*cod, *ncodemax, *Tage; | int firstpass=0, lastpass=4,*cod, *ncodemax, *Tage,*cens; |
| double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint; | double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint; |
| double **pmmij, ***probs; | double **pmmij, ***probs; |
| double *ageexmed,*agecens; | |
| double dateintmean=0; | double dateintmean=0; |
| double *weight; | 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 302 double ftolhess; /* Tolerance for comput | Line 384 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) | |
| the name of the file (name), its extension only (ext) and its first part of the name (finame) | |
| */ | |
| char *ss; /* pointer */ | char *ss; /* pointer */ |
| int l1, l2; /* length counters */ | int l1, l2; /* length counters */ |
| Line 333 static int split( char *path, char *dirc | Line 418 static int split( char *path, char *dirc |
| #endif | #endif |
| */ | */ |
| ss = strrchr( name, '.' ); /* find last / */ | ss = strrchr( name, '.' ); /* find last / */ |
| ss++; | if (ss >0){ |
| strcpy(ext,ss); /* save extension */ | ss++; |
| l1= strlen( name); | strcpy(ext,ss); /* save extension */ |
| l2= strlen(ss)+1; | l1= strlen( name); |
| strncpy( finame, name, l1-l2); | l2= strlen(ss)+1; |
| finame[l1-l2]= 0; | strncpy( finame, name, l1-l2); |
| finame[l1-l2]= 0; | |
| } | |
| return( 0 ); /* we're done */ | return( 0 ); /* we're done */ |
| } | } |
| Line 371 int nbocc(char *s, char occ) | Line 458 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 553 void free_ma3x(double ***m, long nrl, lo | Line 640 void free_ma3x(double ***m, long nrl, lo |
| free((FREE_ARG)(m+nrl-NR_END)); | free((FREE_ARG)(m+nrl-NR_END)); |
| } | } |
| /*************** function subdirf ***********/ | |
| char *subdirf(char fileres[]) | |
| { | |
| /* Caution optionfilefiname is hidden */ | |
| strcpy(tmpout,optionfilefiname); | |
| strcat(tmpout,"/"); /* Add to the right */ | |
| strcat(tmpout,fileres); | |
| return tmpout; | |
| } | |
| /*************** function subdirf2 ***********/ | |
| char *subdirf2(char fileres[], char *preop) | |
| { | |
| /* Caution optionfilefiname is hidden */ | |
| strcpy(tmpout,optionfilefiname); | |
| strcat(tmpout,"/"); | |
| strcat(tmpout,preop); | |
| strcat(tmpout,fileres); | |
| return tmpout; | |
| } | |
| /*************** function subdirf3 ***********/ | |
| char *subdirf3(char fileres[], char *preop, char *preop2) | |
| { | |
| /* Caution optionfilefiname is hidden */ | |
| strcpy(tmpout,optionfilefiname); | |
| strcat(tmpout,"/"); | |
| strcat(tmpout,preop); | |
| strcat(tmpout,preop2); | |
| strcat(tmpout,fileres); | |
| return tmpout; | |
| } | |
| /***************** f1dim *************************/ | /***************** f1dim *************************/ |
| extern int ncom; | extern int ncom; |
| extern double *pcom,*xicom; | extern double *pcom,*xicom; |
| Line 766 void powell(double p[], double **xi, int | Line 888 void powell(double p[], double **xi, int |
| last_time=curr_time; | last_time=curr_time; |
| (void) gettimeofday(&curr_time,&tzp); | (void) gettimeofday(&curr_time,&tzp); |
| printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout); | printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout); |
| fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); | /* fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); |
| fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); | fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); |
| for (i=1;i<=n;i++) { | */ |
| for (i=1;i<=n;i++) { | |
| printf(" %d %.12f",i, p[i]); | printf(" %d %.12f",i, p[i]); |
| fprintf(ficlog," %d %.12lf",i, p[i]); | fprintf(ficlog," %d %.12lf",i, p[i]); |
| fprintf(ficrespow," %.12lf", p[i]); | fprintf(ficrespow," %.12lf", p[i]); |
| Line 778 void powell(double p[], double **xi, int | Line 901 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); |
| asctime_r(&tm,strcurr); | strcpy(strcurr,asctime(&tm)); |
| forecast_time=curr_time; | /* asctime_r(&tm,strcurr); */ |
| 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); |
| for(niterf=10;niterf<=30;niterf+=10){ | for(niterf=10;niterf<=30;niterf+=10){ |
| forecast_time.tv_sec=curr_time.tv_sec+(niterf-*iter)*(curr_time.tv_sec-last_time.tv_sec); | forecast_time.tv_sec=curr_time.tv_sec+(niterf-*iter)*(curr_time.tv_sec-last_time.tv_sec); |
| tmf = *localtime(&forecast_time.tv_sec); | tmf = *localtime(&forecast_time.tv_sec); |
| asctime_r(&tmf,strfor); | /* asctime_r(&tmf,strfor); */ |
| /* strcpy(strfor,asctime(&tmf)); */ | strcpy(strfor,asctime(&tmf)); |
| 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 959 double **pmij(double **ps, double *cov, | Line 1083 double **pmij(double **ps, double *cov, |
| int i,j,j1, nc, ii, jj; | int i,j,j1, nc, ii, jj; |
| for(i=1; i<= nlstate; i++){ | for(i=1; i<= nlstate; i++){ |
| for(j=1; j<i;j++){ | for(j=1; j<i;j++){ |
| for (nc=1, s2=0.;nc <=ncovmodel; nc++){ | for (nc=1, s2=0.;nc <=ncovmodel; nc++){ |
| /*s2 += param[i][j][nc]*cov[nc];*/ | /*s2 += param[i][j][nc]*cov[nc];*/ |
| s2 += x[(i-1)*nlstate*ncovmodel+(j-1)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc]; | s2 += x[(i-1)*nlstate*ncovmodel+(j-1)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc]; |
| /*printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2);*/ | /* printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2); */ |
| } | } |
| ps[i][j]=s2; | ps[i][j]=s2; |
| /*printf("s1=%.17e, s2=%.17e\n",s1,s2);*/ | /* printf("s1=%.17e, s2=%.17e\n",s1,s2); */ |
| } | } |
| for(j=i+1; j<=nlstate+ndeath;j++){ | for(j=i+1; j<=nlstate+ndeath;j++){ |
| for (nc=1, s2=0.;nc <=ncovmodel; nc++){ | for (nc=1, s2=0.;nc <=ncovmodel; nc++){ |
| s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc]; | s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc]; |
| /*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/ | /* printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2); */ |
| } | |
| ps[i][j]=s2; | |
| } | } |
| ps[i][j]=s2; | |
| } | } |
| } | |
| /*ps[3][2]=1;*/ | /*ps[3][2]=1;*/ |
| for(i=1; i<= nlstate; i++){ | for(i=1; i<= nlstate; i++){ |
| s1=0; | s1=0; |
| for(j=1; j<i; j++) | for(j=1; j<i; j++) |
| s1+=exp(ps[i][j]); | s1+=exp(ps[i][j]); |
| for(j=i+1; j<=nlstate+ndeath; j++) | for(j=i+1; j<=nlstate+ndeath; j++) |
| s1+=exp(ps[i][j]); | s1+=exp(ps[i][j]); |
| ps[i][i]=1./(s1+1.); | ps[i][i]=1./(s1+1.); |
| for(j=1; j<i; j++) | for(j=1; j<i; j++) |
| ps[i][j]= exp(ps[i][j])*ps[i][i]; | ps[i][j]= exp(ps[i][j])*ps[i][i]; |
| for(j=i+1; j<=nlstate+ndeath; j++) | for(j=i+1; j<=nlstate+ndeath; j++) |
| ps[i][j]= exp(ps[i][j])*ps[i][i]; | ps[i][j]= exp(ps[i][j])*ps[i][i]; |
| /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */ | /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */ |
| } /* end i */ | } /* end i */ |
| for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){ | for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){ |
| for(jj=1; jj<= nlstate+ndeath; jj++){ | for(jj=1; jj<= nlstate+ndeath; jj++){ |
| ps[ii][jj]=0; | ps[ii][jj]=0; |
| ps[ii][ii]=1; | ps[ii][ii]=1; |
| } | |
| } | } |
| } | |
| /* for(ii=1; ii<= nlstate+ndeath; ii++){ | /* for(ii=1; ii<= nlstate+ndeath; ii++){ */ |
| for(jj=1; jj<= nlstate+ndeath; jj++){ | /* for(jj=1; jj<= nlstate+ndeath; jj++){ */ |
| printf("%lf ",ps[ii][jj]); | /* printf("ddd %lf ",ps[ii][jj]); */ |
| } | /* } */ |
| printf("\n "); | /* printf("\n "); */ |
| } | /* } */ |
| printf("\n ");printf("%lf ",cov[2]);*/ | /* printf("\n ");printf("%lf ",cov[2]); */ |
| /* | /* |
| for(i=1; i<= npar; i++) printf("%f ",x[i]); | for(i=1; i<= npar; i++) printf("%f ",x[i]); |
| goto end;*/ | goto end;*/ |
| return ps; | return ps; |
| } | } |
| Line 1133 double func( double *x) | Line 1257 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 1148 double func( double *x) | Line 1272 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 1176 double func( double *x) | Line 1301 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 1209 double func( double *x) | Line 1355 double func( double *x) |
| oldm=newm; | oldm=newm; |
| } /* end mult */ | } /* end mult */ |
| /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ | |
| /* But now since version 0.9 we anticipate for bias and large stepm. | |
| * If stepm is larger than one month (smallest stepm) and if the exact delay | |
| * (in months) between two waves is not a multiple of stepm, we rounded to | |
| * the nearest (and in case of equal distance, to the lowest) interval but now | |
| * we keep into memory the bias bh[mi][i] and also the previous matrix product | |
| * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the | |
| * probability in order to take into account the bias as a fraction of the way | |
| * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies | |
| * -stepm/2 to stepm/2 . | |
| * For stepm=1 the results are the same as for previous versions of Imach. | |
| * For stepm > 1 the results are less biased than in previous versions. | |
| */ | |
| s1=s[mw[mi][i]][i]; | 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 | |
| * is higher than the multiple of stepm and negative otherwise. | |
| */ | |
| lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ | lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ |
| /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ | |
| /*lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-+bh)*out[s1][s2])); */ /* exponential interpolation */ | |
| /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ | |
| /*if(lli ==000.0)*/ | |
| /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ | |
| ipmx +=1; | ipmx +=1; |
| sw += weight[i]; | sw += weight[i]; |
| ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; | ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
| Line 1260 double func( double *x) | Line 1385 double func( double *x) |
| oldm=newm; | oldm=newm; |
| } /* end mult */ | } /* end mult */ |
| /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ | |
| /* But now since version 0.9 we anticipate for bias and large stepm. | |
| * If stepm is larger than one month (smallest stepm) and if the exact delay | |
| * (in months) between two waves is not a multiple of stepm, we rounded to | |
| * the nearest (and in case of equal distance, to the lowest) interval but now | |
| * we keep into memory the bias bh[mi][i] and also the previous matrix product | |
| * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the | |
| * probability in order to take into account the bias as a fraction of the way | |
| * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies | |
| * -stepm/2 to stepm/2 . | |
| * For stepm=1 the results are the same as for previous versions of Imach. | |
| * For stepm > 1 the results are less biased than in previous versions. | |
| */ | |
| s1=s[mw[mi][i]][i]; | 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 | |
| * is higher than the multiple of stepm and negative otherwise. | |
| */ | |
| /* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); */ /* linear interpolation */ | |
| lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ | lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ |
| /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ | |
| /*if(lli ==000.0)*/ | |
| /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ | |
| ipmx +=1; | ipmx +=1; |
| sw += weight[i]; | sw += weight[i]; |
| ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; | ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
| Line 1449 double funcone( double *x) | Line 1554 double funcone( double *x) |
| return -l; | return -l; |
| } | } |
| char *subdirf(char fileres[]) | |
| { | |
| /* Caution optionfilefiname is hidden */ | |
| strcpy(tmpout,optionfilefiname); | |
| strcat(tmpout,"/"); /* Add to the right */ | |
| strcat(tmpout,fileres); | |
| return tmpout; | |
| } | |
| char *subdirf2(char fileres[], char *preop) | |
| { | |
| strcpy(tmpout,optionfilefiname); | |
| strcat(tmpout,"/"); | |
| strcat(tmpout,preop); | |
| strcat(tmpout,fileres); | |
| return tmpout; | |
| } | |
| char *subdirf3(char fileres[], char *preop, char *preop2) | |
| { | |
| strcpy(tmpout,optionfilefiname); | |
| strcat(tmpout,"/"); | |
| strcat(tmpout,preop); | |
| strcat(tmpout,preop2); | |
| strcat(tmpout,fileres); | |
| return tmpout; | |
| } | |
| /*************** function likelione ***********/ | |
| void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*funcone)(double [])) | void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*funcone)(double [])) |
| { | { |
| /* This routine should help understanding what is done with | /* This routine should help understanding what is done with |
| Line 1520 void mlikeli(FILE *ficres,double p[], in | Line 1598 void mlikeli(FILE *ficres,double p[], in |
| double **xi; | double **xi; |
| double fret; | double fret; |
| double fretone; /* Only one call to likelihood */ | double fretone; /* Only one call to likelihood */ |
| char filerespow[FILENAMELENGTH]; | /* char filerespow[FILENAMELENGTH];*/ |
| xi=matrix(1,npar,1,npar); | xi=matrix(1,npar,1,npar); |
| for (i=1;i<=npar;i++) | for (i=1;i<=npar;i++) |
| for (j=1;j<=npar;j++) | for (j=1;j<=npar;j++) |
| Line 1555 void hesscov(double **matcov, double p[] | Line 1633 void hesscov(double **matcov, double p[] |
| int i, j,jk; | int i, j,jk; |
| int *indx; | int *indx; |
| double hessii(double p[], double delta, int theta, double delti[]); | double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar); |
| double hessij(double p[], double delti[], int i, int j); | double hessij(double p[], double delti[], int i, int j,double (*func)(double []),int npar); |
| void lubksb(double **a, int npar, int *indx, double b[]) ; | void lubksb(double **a, int npar, int *indx, double b[]) ; |
| void ludcmp(double **a, int npar, int *indx, double *d) ; | void ludcmp(double **a, int npar, int *indx, double *d) ; |
| double gompertz(double p[]); | |
| hess=matrix(1,npar,1,npar); | hess=matrix(1,npar,1,npar); |
| printf("\nCalculation of the hessian matrix. Wait...\n"); | printf("\nCalculation of the hessian matrix. Wait...\n"); |
| Line 1567 void hesscov(double **matcov, double p[] | Line 1645 void hesscov(double **matcov, double p[] |
| for (i=1;i<=npar;i++){ | for (i=1;i<=npar;i++){ |
| printf("%d",i);fflush(stdout); | printf("%d",i);fflush(stdout); |
| fprintf(ficlog,"%d",i);fflush(ficlog); | fprintf(ficlog,"%d",i);fflush(ficlog); |
| hess[i][i]=hessii(p,ftolhess,i,delti); | |
| /*printf(" %f ",p[i]);*/ | hess[i][i]=hessii(p,ftolhess,i,delti,func,npar); |
| /*printf(" %lf ",hess[i][i]);*/ | |
| /* printf(" %f ",p[i]); | |
| printf(" %lf %lf %lf",hess[i][i],ftolhess,delti[i]);*/ | |
| } | } |
| for (i=1;i<=npar;i++) { | for (i=1;i<=npar;i++) { |
| Line 1577 void hesscov(double **matcov, double p[] | Line 1657 void hesscov(double **matcov, double p[] |
| if (j>i) { | if (j>i) { |
| printf(".%d%d",i,j);fflush(stdout); | printf(".%d%d",i,j);fflush(stdout); |
| fprintf(ficlog,".%d%d",i,j);fflush(ficlog); | fprintf(ficlog,".%d%d",i,j);fflush(ficlog); |
| hess[i][j]=hessij(p,delti,i,j); | hess[i][j]=hessij(p,delti,i,j,func,npar); |
| hess[j][i]=hess[i][j]; | hess[j][i]=hess[i][j]; |
| /*printf(" %lf ",hess[i][j]);*/ | /*printf(" %lf ",hess[i][j]);*/ |
| } | } |
| Line 1648 void hesscov(double **matcov, double p[] | Line 1729 void hesscov(double **matcov, double p[] |
| } | } |
| /*************** hessian matrix ****************/ | /*************** hessian matrix ****************/ |
| double hessii( double x[], double delta, int theta, double delti[]) | double hessii(double x[], double delta, int theta, double delti[], double (*func)(double []), int npar) |
| { | { |
| int i; | int i; |
| int l=1, lmax=20; | int l=1, lmax=20; |
| double k1,k2; | double k1,k2; |
| double p2[NPARMAX+1]; | double p2[NPARMAX+1]; |
| double res; | double res; |
| double delt, delts, nkhi=10.,nkhif=1., khi=1.e-4; | double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; |
| double fx; | double fx; |
| int k=0,kmax=10; | int k=0,kmax=10; |
| double l1; | double l1; |
| Line 1695 double hessii( double x[], double delta, | Line 1776 double hessii( double x[], double delta, |
| } | } |
| double hessij( double x[], double delti[], int thetai,int thetaj) | double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) |
| { | { |
| int i; | int i; |
| int l=1, l1, lmax=20; | int l=1, l1, lmax=20; |
| Line 1805 void lubksb(double **a, int n, int *indx | Line 1886 void lubksb(double **a, int n, int *indx |
| } | } |
| /************ 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 1825 void freqsummary(char fileres[], int ia | Line 1906 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 1838 void freqsummary(char fileres[], int ia | Line 1919 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 1878 void freqsummary(char fileres[], int ia | Line 1959 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);*/ |
| fprintf(ficresp, "#Local time at start: %s", strstart); | |
| 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 1964 void freqsummary(char fileres[], int ia | Line 2045 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 2073 void concatwav(int wav[], int **dh, int | Line 2154 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 2113 void concatwav(int wav[], int **dh, int | Line 2194 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) jmax=j; |
| Line 2127 void concatwav(int wav[], int **dh, int | Line 2208 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) jmax=j; |
| else if (j <= jmin)jmin=j; | else if (j <= jmin)jmin=j; |
| Line 2220 void tricode(int *Tvar, int **nbcode, in | Line 2302 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 2238 void tricode(int *Tvar, int **nbcode, in | Line 2320 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 ij, int estepm,double delti[],double **matcov,char strstart[] ) |
| { | { |
| /* Health expectancies */ | /* Health expectancies */ |
| Line 2256 void evsij(char fileres[], double ***eij | Line 2338 void evsij(char fileres[], double ***eij |
| 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,"# Local time at start: %s", strstart); | |
| fprintf(ficreseij,"# Health expectancies\n"); | fprintf(ficreseij,"# Health expectancies\n"); |
| fprintf(ficreseij,"# Age"); | fprintf(ficreseij,"# Age"); |
| for(i=1; i<=nlstate;i++) | for(i=1; i<=nlstate;i++) |
| Line 2311 void evsij(char fileres[], double ***eij | Line 2394 void evsij(char fileres[], double ***eij |
| 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 */ |
| for(theta=1; theta <=npar; theta++){ | for(theta=1; theta <=npar; theta++){ |
| for(i=1; i<=npar; i++){ | for(i=1; i<=npar; i++){ |
| Line 2410 void evsij(char fileres[], double ***eij | Line 2493 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 2461 void varevsij(char optionfilefiname[], d | Line 2544 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); |
| fprintf(ficresprobmorprev, "#Local time at start: %s", strstart); | |
| 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 2471 void varevsij(char optionfilefiname[], d | Line 2556 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); |
| fprintf(ficresvij, "#Local time at start: %s", strstart); | |
| 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 the stable prevalence in health states i\n"); |
| fprintf(ficresvij,"# Age"); | fprintf(ficresvij,"# Age"); |
| for(i=1; i<=nlstate;i++) | for(i=1; i<=nlstate;i++) |
| Line 2708 void varevsij(char optionfilefiname[], d | Line 2794 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 2721 void varprevlim(char fileres[], double * | Line 2807 void varprevlim(char fileres[], double * |
| double **gradg, **trgradg; | double **gradg, **trgradg; |
| double age,agelim; | double age,agelim; |
| int theta; | int theta; |
| fprintf(ficresvpl, "#Local time at start: %s", strstart); | |
| fprintf(ficresvpl,"# Standard deviation of stable prevalences \n"); | fprintf(ficresvpl,"# Standard deviation of stable prevalences \n"); |
| fprintf(ficresvpl,"# Age"); | fprintf(ficresvpl,"# Age"); |
| for(i=1; i<=nlstate;i++) | for(i=1; i<=nlstate;i++) |
| Line 2791 void varprevlim(char fileres[], double * | Line 2877 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 2836 void varprob(char optionfilefiname[], do | Line 2922 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); |
| fprintf(ficresprob, "#Local time at start: %s", strstart); | |
| 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"); |
| fprintf(ficresprobcov, "#Local time at start: %s", strstart); | |
| fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); | fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); |
| fprintf(ficresprobcov,"# Age"); | fprintf(ficresprobcov,"# Age"); |
| fprintf(ficresprobcor, "#Local time at start: %s", strstart); | |
| fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); | fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); |
| fprintf(ficresprobcov,"# Age"); | fprintf(ficresprobcov,"# Age"); |
| Line 2865 void varprob(char optionfilefiname[], do | Line 2953 void varprob(char optionfilefiname[], do |
| fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n"); | fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n"); |
| fprintf(fichtm,"\n"); | fprintf(fichtm,"\n"); |
| fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Computing matrix of variance-covariance of step probabilities</a></h4></li>\n",optionfilehtmcov); | fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of pairs of step probabilities (drawings)</a></h4></li>\n",optionfilehtmcov); |
| fprintf(fichtmcov,"\n<h4>Computing matrix of variance-covariance of step probabilities</h4>\n\ | fprintf(fichtmcov,"\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n\ |
| file %s<br>\n",optionfilehtmcov); | file %s<br>\n",optionfilehtmcov); |
| fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated\ | fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated\ |
| and drawn. It helps understanding how is the covariance between two incidences.\ | and drawn. It helps understanding how is the covariance between two incidences.\ |
| Line 2899 To be simple, these graphs help to under | Line 2987 To be simple, these graphs help to under |
| fprintf(ficgp, "**********\n#\n"); | fprintf(ficgp, "**********\n#\n"); |
| fprintf(fichtm, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable "); | fprintf(fichtmcov, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable "); |
| for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); | for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); |
| fprintf(fichtm, "**********\n<hr size=\"2\" color=\"#EC5E5E\">"); | fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">"); |
| fprintf(ficresprobcor, "\n#********** Variable "); | fprintf(ficresprobcor, "\n#********** Variable "); |
| for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); | for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); |
| Line 3110 void printinghtml(char fileres[], char t | Line 3198 void printinghtml(char fileres[], char t |
| double jprev1, double mprev1,double anprev1, \ | double jprev1, double mprev1,double anprev1, \ |
| double jprev2, double mprev2,double anprev2){ | double jprev2, double mprev2,double anprev2){ |
| int jj1, k1, i1, cpt; | int jj1, k1, i1, cpt; |
| /*char optionfilehtm[FILENAMELENGTH];*/ | |
| /* if((fichtm=fopen(optionfilehtm,"a"))==NULL) { */ | |
| /* printf("Problem with %s \n",optionfilehtm), exit(0); */ | |
| /* fprintf(ficlog,"Problem with %s \n",optionfilehtm), exit(0); */ | |
| /* } */ | |
| fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n \ | fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</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 ", |
| - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n \ | jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p")); |
| - Stable prevalence in each health state: <a href=\"%s\">%s</a> <br>\n \ | fprintf(fichtm,"\ |
| - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ", | |
| stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij")); | |
| fprintf(fichtm,"\ | |
| - Stable prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", | |
| subdirf2(fileres,"pl"),subdirf2(fileres,"pl")); | |
| fprintf(fichtm,"\ | |
| - Life expectancies by age and initial health status (estepm=%2d months): \ | - Life expectancies by age and initial health status (estepm=%2d months): \ |
| <a href=\"%s\">%s</a> <br>\n</li>", \ | <a href=\"%s\">%s</a> <br>\n</li>", |
| jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p"),\ | |
| stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij"),\ | |
| subdirf2(fileres,"pl"),subdirf2(fileres,"pl"),\ | |
| 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>"); |
| Line 3158 fprintf(fichtm," \n<ul><li><b>Graphs</b> | Line 3244 fprintf(fichtm," \n<ul><li><b>Graphs</b> |
| fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): %s%d%d.png <br> \ | fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): %s%d%d.png <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); |
| } | } |
| fprintf(fichtm,"\n<br>- Total life expectancy by age and \ | |
| health expectancies in states (1) and (2): %s%d.png<br>\ | |
| <img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1); | |
| } /* end i1 */ | } /* end i1 */ |
| }/* End k1 */ | }/* End k1 */ |
| fprintf(fichtm,"</ul>"); | fprintf(fichtm,"</ul>"); |
| fprintf(fichtm,"\n<br><li><h4> Result files (second order: variances)</h4>\n\ | fprintf(fichtm,"\ |
| - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n\ | \n<br><li><h4> Result files (second order: variances)</h4>\n\ |
| - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n\ | - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n", rfileres,rfileres); |
| - Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n\ | |
| - Correlation matrix 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", |
| - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"%s\">%s</a><br>\n\ | subdirf2(fileres,"prob"),subdirf2(fileres,"prob")); |
| - Health expectancies with their variances (no covariance): <a href=\"%s\">%s</a> <br>\n\ | fprintf(fichtm,"\ |
| - Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", | |
| subdirf2(fileres,"probcov"),subdirf2(fileres,"probcov")); | |
| fprintf(fichtm,"\ | |
| - Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", | |
| subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor")); | |
| fprintf(fichtm,"\ | |
| - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"%s\">%s</a><br>\n", | |
| estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v")); | |
| fprintf(fichtm,"\ | |
| - Health expectancies with their variances (no covariance): <a href=\"%s\">%s</a> <br>\n", | |
| subdirf2(fileres,"t"),subdirf2(fileres,"t")); | |
| fprintf(fichtm,"\ | |
| - Standard deviation of stable prevalences: <a href=\"%s\">%s</a> <br>\n",\ | - Standard deviation of stable prevalences: <a href=\"%s\">%s</a> <br>\n",\ |
| rfileres,rfileres,\ | |
| subdirf2(fileres,"prob"),subdirf2(fileres,"prob"),\ | |
| subdirf2(fileres,"probcov"),subdirf2(fileres,"probcov"),\ | |
| subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor"),\ | |
| estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"),\ | |
| subdirf2(fileres,"t"),subdirf2(fileres,"t"),\ | |
| subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl")); | subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl")); |
| /* if(popforecast==1) fprintf(fichtm,"\n */ | /* if(popforecast==1) fprintf(fichtm,"\n */ |
| Line 3188 health expectancies in states (1) and (2 | Line 3278 health expectancies in states (1) and (2 |
| /* <br>",fileres,fileres,fileres,fileres); */ | /* <br>",fileres,fileres,fileres,fileres); */ |
| /* else */ | /* else */ |
| /* fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */ | /* fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */ |
| fprintf(fichtm," <ul><li><b>Graphs</b></li><p>"); | fflush(fichtm); |
| fprintf(fichtm," <ul><li><b>Graphs</b></li><p>"); | |
| m=cptcoveff; | m=cptcoveff; |
| if (cptcovn < 1) {m=1;ncodemax[1]=1;} | if (cptcovn < 1) {m=1;ncodemax[1]=1;} |
| Line 3204 fprintf(fichtm," <ul><li><b>Graphs</b></ | Line 3295 fprintf(fichtm," <ul><li><b>Graphs</b></ |
| fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); | fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
| } | } |
| for(cpt=1; cpt<=nlstate;cpt++) { | for(cpt=1; cpt<=nlstate;cpt++) { |
| fprintf(fichtm,"<br>- Observed and period prevalence (with confident\ | fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \ |
| interval) in state (%d): %s%d%d.png <br>\ | prevalence (with 95%% confidence interval) in state (%d): %s%d%d.png <br>\ |
| <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1); | <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1); |
| } | } |
| fprintf(fichtm,"\n<br>- Total life expectancy by age and \ | |
| health expectancies in states (1) and (2): %s%d.png<br>\ | |
| <img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1); | |
| } /* end i1 */ | } /* end i1 */ |
| }/* End k1 */ | }/* End k1 */ |
| fprintf(fichtm,"</ul>"); | fprintf(fichtm,"</ul>"); |
| Line 3741 int fileappend(FILE *fichier, char *opti | Line 3835 int fileappend(FILE *fichier, char *opti |
| fflush(fichier); | fflush(fichier); |
| return (1); | return (1); |
| } | } |
| /**************** function prwizard **********************/ | |
| void prwizard(int ncovmodel, int nlstate, int ndeath, char model[], FILE *ficparo) | void prwizard(int ncovmodel, int nlstate, int ndeath, char model[], FILE *ficparo) |
| { | { |
| /* Wizard to print covariance matrix template */ | |
| char ca[32], cb[32], cc[32]; | char ca[32], cb[32], cc[32]; |
| int i,j, k, l, li, lj, lk, ll, jj, npar, itimes; | int i,j, k, l, li, lj, lk, ll, jj, npar, itimes; |
| int numlinepar; | int numlinepar; |
| Line 3857 void prwizard(int ncovmodel, int nlstate | Line 3956 void prwizard(int ncovmodel, int nlstate |
| } /* end k*/ | } /* end k*/ |
| } /*end j */ | } /*end j */ |
| } /* end i */ | } /* end i */ |
| } | } /* end itimes */ |
| } /* end of prwizard */ | } /* end of prwizard */ |
| /******************* Gompertz Likelihood ******************************/ | |
| double gompertz(double x[]) | |
| { | |
| double A,B,L=0.0,sump=0.,num=0.; | |
| int i,n=0; /* n is the size of the sample */ | |
| for (i=0;i<=imx-1 ; i++) { | |
| sump=sump+weight[i]; | |
| /* sump=sump+1;*/ | |
| num=num+1; | |
| } | |
| /* for (i=0; i<=imx; i++) | |
| if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/ | |
| for (i=1;i<=imx ; i++) | |
| { | |
| if (cens[i]==1 & wav[i]>1) | |
| A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); | |
| if (cens[i]==0 & wav[i]>1) | |
| A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))) | |
| +log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM); | |
| if (wav[i]>1 & agecens[i]>15) { | |
| L=L+A*weight[i]; | |
| /* printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/ | |
| } | |
| } | |
| /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/ | |
| return -2*L*num/sump; | |
| } | |
| /******************* Printing html file ***********/ | |
| void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \ | |
| int lastpass, int stepm, int weightopt, char model[],\ | |
| int imx, double p[],double **matcov,double agemortsup){ | |
| int i,k; | |
| fprintf(fichtm,"<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); | |
| for (i=1;i<=2;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,"</ul>"); | |
| fprintf(fichtm,"<ul><li><h4>Life table</h4>\n <br>"); | |
| fprintf(fichtm,"\nAge lx qx dx Lx Tx e(x)<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); | |
| } | |
| /******************* Gnuplot file **************/ | |
| void printinggnuplotmort(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ | |
| char dirfileres[132],optfileres[132]; | |
| int m,cpt,k1,i,k,j,jk,k2,k3,ij,l; | |
| int ng; | |
| /*#ifdef windows */ | |
| fprintf(ficgp,"cd \"%s\" \n",pathc); | |
| /*#endif */ | |
| strcpy(dirfileres,optionfilefiname); | |
| strcpy(optfileres,"vpl"); | |
| fprintf(ficgp,"set out \"graphmort.png\"\n "); | |
| fprintf(ficgp,"set xlabel \"Age\"\n set ylabel \"Force of mortality (per year)\" \n "); | |
| fprintf(ficgp, "set ter png small\n set log y\n"); | |
| fprintf(ficgp, "set size 0.65,0.65\n"); | |
| fprintf(ficgp,"plot [%d:100] %lf*exp(%lf*(x-%d))",agegomp,p[1],p[2],agegomp); | |
| } | |
| /***********************************************/ | /***********************************************/ |
| Line 3870 int main(int argc, char *argv[]) | Line 4052 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 jj, 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 NDIM=2; | |
| char ca[32], cb[32], cc[32]; | |
| /* FILE *fichtm; *//* Html File */ | /* FILE *fichtm; *//* Html File */ |
| /* FILE *ficgp;*/ /*Gnuplot File */ | /* FILE *ficgp;*/ /*Gnuplot File */ |
| double agedeb, agefin,hf; | double agedeb, agefin,hf; |
| Line 3886 int main(int argc, char *argv[]) | Line 4072 int main(int argc, char *argv[]) |
| int *indx; | int *indx; |
| 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]; | char pathr[MAXLINE], pathimach[MAXLINE]; |
| 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 3896 int main(int argc, char *argv[]) | Line 4082 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 3913 int main(int argc, char *argv[]) | Line 4101 int main(int argc, char *argv[]) |
| double *epj, vepp; | double *epj, vepp; |
| double kk1, kk2; | double kk1, kk2; |
| double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000; | double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000; |
| double **ximort; | |
| char *alph[]={"a","a","b","c","d","e"}, str[4]; | char *alph[]={"a","a","b","c","d","e"}, str[4]; |
| int *dcwave; | |
| char z[1]="c", occ; | char z[1]="c", occ; |
| Line 3926 int main(int argc, char *argv[]) | Line 4114 int main(int argc, char *argv[]) |
| long total_usecs; | long total_usecs; |
| /* setlocale (LC_ALL, ""); */ | |
| /* bindtextdomain (PACKAGE, LOCALEDIR); */ | |
| /* textdomain (PACKAGE); */ | |
| /* setlocale (LC_CTYPE, ""); */ | |
| /* setlocale (LC_MESSAGES, ""); */ | |
| /* gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ | /* gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ |
| (void) gettimeofday(&start_time,&tzp); | (void) gettimeofday(&start_time,&tzp); |
| curr_time=start_time; | curr_time=start_time; |
| Line 3966 int main(int argc, char *argv[]) | Line 4160 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],pathimach,optionfile,optionfilext,optionfilefiname); | |
| /* strcpy(pathimach,argv[0]); */ | |
| split(pathtot,path,optionfile,optionfilext,optionfilefiname); | split(pathtot,path,optionfile,optionfilext,optionfilefiname); |
| printf("pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); | printf("pathimach=%s, pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname); |
| chdir(path); | chdir(path); |
| strcpy(command,"mkdir "); | strcpy(command,"mkdir "); |
| strcat(command,optionfilefiname); | strcat(command,optionfilefiname); |
| Line 3992 int main(int argc, char *argv[]) | Line 4188 int main(int argc, char *argv[]) |
| } | } |
| fprintf(ficlog,"Log filename:%s\n",filelog); | fprintf(ficlog,"Log filename:%s\n",filelog); |
| fprintf(ficlog,"\n%s\n%s",version,fullversion); | fprintf(ficlog,"\n%s\n%s",version,fullversion); |
| fprintf(ficlog,"\nEnter the parameter file name: "); | fprintf(ficlog,"\nEnter the parameter file name: \n"); |
| fprintf(ficlog,"pathtot=%s\n\ | fprintf(ficlog,"pathimach=%s\npathtot=%s\n\ |
| path=%s \n\ | path=%s \n\ |
| optionfile=%s\n\ | optionfile=%s\n\ |
| optionfilext=%s\n\ | optionfilext=%s\n\ |
| optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); | optionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname); |
| printf("Localtime (at start):%s",strstart); | printf("Local time (at start):%s",strstart); |
| fprintf(ficlog,"Localtime (at start): %s",strstart); | fprintf(ficlog,"Local time (at start): %s",strstart); |
| fflush(ficlog); | fflush(ficlog); |
| /* (void) gettimeofday(&curr_time,&tzp); */ | /* (void) gettimeofday(&curr_time,&tzp); */ |
| /* printf("Elapsed time %d\n", asc_diff_time(curr_time.tv_sec-start_time.tv_sec,tmpout)); */ | /* printf("Elapsed time %d\n", asc_diff_time(curr_time.tv_sec-start_time.tv_sec,tmpout)); */ |
| Line 4065 int main(int argc, char *argv[]) | Line 4261 int main(int argc, char *argv[]) |
| ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ | ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ |
| nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ | nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ |
| npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ | |
| delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); | |
| delti=delti3[1][1]; | |
| /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ | |
| if(mle==-1){ /* Print a wizard for help writing covariance matrix */ | if(mle==-1){ /* Print a wizard for help writing covariance matrix */ |
| prwizard(ncovmodel, nlstate, ndeath, model, ficparo); | prwizard(ncovmodel, nlstate, ndeath, model, ficparo); |
| printf(" You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso); | printf(" You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso); |
| fprintf(ficlog," You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso); | fprintf(ficlog," You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso); |
| free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); | |
| fclose (ficparo); | fclose (ficparo); |
| fclose (ficlog); | fclose (ficlog); |
| exit(0); | exit(0); |
| } | } |
| /* Read guess parameters */ | else if(mle==-3) { |
| /* Reads comments: lines beginning with '#' */ | prwizard(ncovmodel, nlstate, ndeath, model, ficparo); |
| while((c=getc(ficpar))=='#' && c!= EOF){ | printf(" You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso); |
| ungetc(c,ficpar); | fprintf(ficlog," You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso); |
| fgets(line, MAXLINE, ficpar); | param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); |
| numlinepar++; | matcov=matrix(1,npar,1,npar); |
| puts(line); | |
| fputs(line,ficparo); | |
| fputs(line,ficlog); | |
| } | } |
| ungetc(c,ficpar); | else{ |
| /* Read guess parameters */ | |
| /* Reads comments: lines beginning with '#' */ | |
| while((c=getc(ficpar))=='#' && c!= EOF){ | |
| ungetc(c,ficpar); | |
| fgets(line, MAXLINE, ficpar); | |
| numlinepar++; | |
| puts(line); | |
| fputs(line,ficparo); | |
| fputs(line,ficlog); | |
| } | |
| ungetc(c,ficpar); | |
| param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); | |
| for(i=1; i <=nlstate; i++){ | |
| j=0; | |
| for(jj=1; jj <=nlstate+ndeath; jj++){ | |
| if(jj==i) continue; | |
| j++; | |
| fscanf(ficpar,"%1d%1d",&i1,&j1); | |
| if ((i1 != i) && (j1 != j)){ | |
| printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); | |
| exit(1); | |
| } | |
| fprintf(ficparo,"%1d%1d",i1,j1); | |
| if(mle==1) | |
| printf("%1d%1d",i,j); | |
| fprintf(ficlog,"%1d%1d",i,j); | |
| for(k=1; k<=ncovmodel;k++){ | |
| fscanf(ficpar," %lf",¶m[i][j][k]); | |
| if(mle==1){ | |
| printf(" %lf",param[i][j][k]); | |
| fprintf(ficlog," %lf",param[i][j][k]); | |
| } | |
| else | |
| fprintf(ficlog," %lf",param[i][j][k]); | |
| fprintf(ficparo," %lf",param[i][j][k]); | |
| } | |
| fscanf(ficpar,"\n"); | |
| numlinepar++; | |
| if(mle==1) | |
| printf("\n"); | |
| fprintf(ficlog,"\n"); | |
| fprintf(ficparo,"\n"); | |
| } | |
| } | |
| fflush(ficlog); | |
| param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); | |
| for(i=1; i <=nlstate; i++){ | p=param[1][1]; |
| j=0; | |
| for(jj=1; jj <=nlstate+ndeath; jj++){ | /* Reads comments: lines beginning with '#' */ |
| if(jj==i) continue; | while((c=getc(ficpar))=='#' && c!= EOF){ |
| j++; | ungetc(c,ficpar); |
| fscanf(ficpar,"%1d%1d",&i1,&j1); | fgets(line, MAXLINE, ficpar); |
| if ((i1 != i) && (j1 != j)){ | numlinepar++; |
| printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); | puts(line); |
| exit(1); | fputs(line,ficparo); |
| fputs(line,ficlog); | |
| } | |
| ungetc(c,ficpar); | |
| for(i=1; i <=nlstate; i++){ | |
| for(j=1; j <=nlstate+ndeath-1; j++){ | |
| fscanf(ficpar,"%1d%1d",&i1,&j1); | |
| if ((i1-i)*(j1-j)!=0){ | |
| printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); | |
| exit(1); | |
| } | |
| printf("%1d%1d",i,j); | |
| fprintf(ficparo,"%1d%1d",i1,j1); | |
| fprintf(ficlog,"%1d%1d",i1,j1); | |
| for(k=1; k<=ncovmodel;k++){ | |
| fscanf(ficpar,"%le",&delti3[i][j][k]); | |
| printf(" %le",delti3[i][j][k]); | |
| fprintf(ficparo," %le",delti3[i][j][k]); | |
| fprintf(ficlog," %le",delti3[i][j][k]); | |
| } | |
| fscanf(ficpar,"\n"); | |
| numlinepar++; | |
| printf("\n"); | |
| fprintf(ficparo,"\n"); | |
| fprintf(ficlog,"\n"); | |
| } | } |
| fprintf(ficparo,"%1d%1d",i1,j1); | } |
| fflush(ficlog); | |
| delti=delti3[1][1]; | |
| /* free_ma3x(delti3,1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */ /* Hasn't to to freed here otherwise delti is no more allocated */ | |
| /* Reads comments: lines beginning with '#' */ | |
| while((c=getc(ficpar))=='#' && c!= EOF){ | |
| ungetc(c,ficpar); | |
| fgets(line, MAXLINE, ficpar); | |
| numlinepar++; | |
| puts(line); | |
| fputs(line,ficparo); | |
| fputs(line,ficlog); | |
| } | |
| ungetc(c,ficpar); | |
| matcov=matrix(1,npar,1,npar); | |
| for(i=1; i <=npar; i++){ | |
| fscanf(ficpar,"%s",&str); | |
| if(mle==1) | if(mle==1) |
| printf("%1d%1d",i,j); | printf("%s",str); |
| fprintf(ficlog,"%1d%1d",i,j); | fprintf(ficlog,"%s",str); |
| for(k=1; k<=ncovmodel;k++){ | fprintf(ficparo,"%s",str); |
| fscanf(ficpar," %lf",¶m[i][j][k]); | for(j=1; j <=i; j++){ |
| fscanf(ficpar," %le",&matcov[i][j]); | |
| if(mle==1){ | if(mle==1){ |
| printf(" %lf",param[i][j][k]); | printf(" %.5le",matcov[i][j]); |
| fprintf(ficlog," %lf",param[i][j][k]); | |
| } | } |
| else | fprintf(ficlog," %.5le",matcov[i][j]); |
| fprintf(ficlog," %lf",param[i][j][k]); | fprintf(ficparo," %.5le",matcov[i][j]); |
| fprintf(ficparo," %lf",param[i][j][k]); | |
| } | } |
| fscanf(ficpar,"\n"); | fscanf(ficpar,"\n"); |
| numlinepar++; | numlinepar++; |
| Line 4118 int main(int argc, char *argv[]) | Line 4407 int main(int argc, char *argv[]) |
| fprintf(ficlog,"\n"); | fprintf(ficlog,"\n"); |
| fprintf(ficparo,"\n"); | fprintf(ficparo,"\n"); |
| } | } |
| } | for(i=1; i <=npar; i++) |
| fflush(ficlog); | for(j=i+1;j<=npar;j++) |
| matcov[i][j]=matcov[j][i]; | |
| npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ | |
| p=param[1][1]; | |
| /* Reads comments: lines beginning with '#' */ | |
| while((c=getc(ficpar))=='#' && c!= EOF){ | |
| ungetc(c,ficpar); | |
| fgets(line, MAXLINE, ficpar); | |
| numlinepar++; | |
| puts(line); | |
| fputs(line,ficparo); | |
| fputs(line,ficlog); | |
| } | |
| ungetc(c,ficpar); | |
| delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); | |
| /* delti=vector(1,npar); *//* Scale of each paramater (output from hesscov) */ | |
| for(i=1; i <=nlstate; i++){ | |
| for(j=1; j <=nlstate+ndeath-1; j++){ | |
| fscanf(ficpar,"%1d%1d",&i1,&j1); | |
| if ((i1-i)*(j1-j)!=0){ | |
| printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); | |
| exit(1); | |
| } | |
| printf("%1d%1d",i,j); | |
| fprintf(ficparo,"%1d%1d",i1,j1); | |
| fprintf(ficlog,"%1d%1d",i1,j1); | |
| for(k=1; k<=ncovmodel;k++){ | |
| fscanf(ficpar,"%le",&delti3[i][j][k]); | |
| printf(" %le",delti3[i][j][k]); | |
| fprintf(ficparo," %le",delti3[i][j][k]); | |
| fprintf(ficlog," %le",delti3[i][j][k]); | |
| } | |
| fscanf(ficpar,"\n"); | |
| numlinepar++; | |
| printf("\n"); | |
| fprintf(ficparo,"\n"); | |
| fprintf(ficlog,"\n"); | |
| } | |
| } | |
| fflush(ficlog); | |
| delti=delti3[1][1]; | |
| /* free_ma3x(delti3,1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */ /* Hasn't to to freed here otherwise delti is no more allocated */ | |
| /* Reads comments: lines beginning with '#' */ | |
| while((c=getc(ficpar))=='#' && c!= EOF){ | |
| ungetc(c,ficpar); | |
| fgets(line, MAXLINE, ficpar); | |
| numlinepar++; | |
| puts(line); | |
| fputs(line,ficparo); | |
| fputs(line,ficlog); | |
| } | |
| ungetc(c,ficpar); | |
| matcov=matrix(1,npar,1,npar); | |
| for(i=1; i <=npar; i++){ | |
| fscanf(ficpar,"%s",&str); | |
| if(mle==1) | |
| printf("%s",str); | |
| fprintf(ficlog,"%s",str); | |
| fprintf(ficparo,"%s",str); | |
| for(j=1; j <=i; j++){ | |
| fscanf(ficpar," %le",&matcov[i][j]); | |
| if(mle==1){ | |
| printf(" %.5le",matcov[i][j]); | |
| } | |
| fprintf(ficlog," %.5le",matcov[i][j]); | |
| fprintf(ficparo," %.5le",matcov[i][j]); | |
| } | |
| fscanf(ficpar,"\n"); | |
| numlinepar++; | |
| if(mle==1) | if(mle==1) |
| printf("\n"); | printf("\n"); |
| fprintf(ficlog,"\n"); | fprintf(ficlog,"\n"); |
| fprintf(ficparo,"\n"); | |
| } | |
| for(i=1; i <=npar; i++) | |
| for(j=i+1;j<=npar;j++) | |
| matcov[i][j]=matcov[j][i]; | |
| if(mle==1) | |
| printf("\n"); | |
| fprintf(ficlog,"\n"); | |
| fflush(ficlog); | |
| /*-------- Rewriting paramater file ----------*/ | |
| strcpy(rfileres,"r"); /* "Rparameterfile */ | |
| strcat(rfileres,optionfilefiname); /* Parameter file first name*/ | |
| strcat(rfileres,"."); /* */ | |
| strcat(rfileres,optionfilext); /* Other files have txt extension */ | |
| if((ficres =fopen(rfileres,"w"))==NULL) { | |
| printf("Problem writing new parameter file: %s\n", fileres);goto end; | |
| fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end; | |
| } | |
| fprintf(ficres,"#%s\n",version); | |
| fflush(ficlog); | |
| /*-------- Rewriting parameter file ----------*/ | |
| strcpy(rfileres,"r"); /* "Rparameterfile */ | |
| strcat(rfileres,optionfilefiname); /* Parameter file first name*/ | |
| strcat(rfileres,"."); /* */ | |
| strcat(rfileres,optionfilext); /* Other files have txt extension */ | |
| if((ficres =fopen(rfileres,"w"))==NULL) { | |
| printf("Problem writing new parameter file: %s\n", fileres);goto end; | |
| fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end; | |
| } | |
| fprintf(ficres,"#%s\n",version); | |
| } /* End of mle != -3 */ | |
| /*-------- data file ----------*/ | /*-------- data file ----------*/ |
| if((fic=fopen(datafile,"r"))==NULL) { | if((fic=fopen(datafile,"r"))==NULL) { |
| printf("Problem with datafile: %s\n", datafile);goto end; | printf("Problem with datafile: %s\n", datafile);goto end; |
| Line 4249 int main(int argc, char *argv[]) | Line 4456 int main(int argc, char *argv[]) |
| i=1; | i=1; |
| while (fgets(line, MAXLINE, fic) != NULL) { | while (fgets(line, MAXLINE, fic) != NULL) { |
| if ((i >= firstobs) && (i <=lastobs)) { | if ((i >= firstobs) && (i <=lastobs)) { |
| for(j=0; line[j] != '\n';j++){ /* Untabifies line */ | |
| if(line[j] == '\t') | |
| line[j] = ' '; | |
| } | |
| for (j=maxwav;j>=1;j--){ | for (j=maxwav;j>=1;j--){ |
| cutv(stra, strb,line,' '); s[j][i]=atoi(strb); | cutv(stra, strb,line,' '); s[j][i]=atoi(strb); |
| strcpy(line,stra); | strcpy(line,stra); |
| Line 4412 int main(int argc, char *argv[]) | Line 4622 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 4472 int main(int argc, char *argv[]) | Line 4682 int main(int argc, char *argv[]) |
| }*/ | }*/ |
| printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); | printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); |
| fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); | fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); |
| agegomp=(int)agemin; | |
| free_vector(severity,1,maxwav); | free_vector(severity,1,maxwav); |
| free_imatrix(outcome,1,maxwav+1,1,n); | free_imatrix(outcome,1,maxwav+1,1,n); |
| free_vector(moisnais,1,n); | free_vector(moisnais,1,n); |
| Line 4528 int main(int argc, char *argv[]) | Line 4740 int main(int argc, char *argv[]) |
| /*------------ gnuplot -------------*/ | /*------------ gnuplot -------------*/ |
| strcpy(optionfilegnuplot,optionfilefiname); | strcpy(optionfilegnuplot,optionfilefiname); |
| if(mle==-3) | |
| strcat(optionfilegnuplot,"-mort"); | |
| strcat(optionfilegnuplot,".gp"); | strcat(optionfilegnuplot,".gp"); |
| if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { | if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { |
| printf("Problem with file %s",optionfilegnuplot); | printf("Problem with file %s",optionfilegnuplot); |
| } | } |
| Line 4541 int main(int argc, char *argv[]) | Line 4756 int main(int argc, char *argv[]) |
| /*--------- index.htm --------*/ | /*--------- index.htm --------*/ |
| strcpy(optionfilehtm,optionfilefiname); /* Main html file */ | strcpy(optionfilehtm,optionfilefiname); /* Main html file */ |
| if(mle==-3) | |
| strcat(optionfilehtm,"-mort"); | |
| strcat(optionfilehtm,".htm"); | strcat(optionfilehtm,".htm"); |
| if((fichtm=fopen(optionfilehtm,"w"))==NULL) { | if((fichtm=fopen(optionfilehtm,"w"))==NULL) { |
| printf("Problem with %s \n",optionfilehtm), exit(0); | printf("Problem with %s \n",optionfilehtm), exit(0); |
| Line 4576 Title=%s <br>Datafile=%s Firstpass=%d La | Line 4793 Title=%s <br>Datafile=%s Firstpass=%d La |
| strcpy(pathr,path); | strcpy(pathr,path); |
| strcat(pathr,optionfilefiname); | strcat(pathr,optionfilefiname); |
| chdir(optionfilefiname); /* Move to directory named optionfile */ | chdir(optionfilefiname); /* Move to directory named optionfile */ |
| strcpy(lfileres,fileres); | |
| strcat(lfileres,"/"); | |
| strcat(lfileres,optionfilefiname); | |
| /* 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 4601 Interval (in months) between two waves: | Line 4815 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*/ |
| likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ | if (mle==-3){ |
| printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); | ximort=matrix(1,NDIM,1,NDIM); |
| for (k=1; k<=npar;k++) | cens=ivector(1,n); |
| printf(" %d %8.5f",k,p[k]); | ageexmed=vector(1,n); |
| printf("\n"); | agecens=vector(1,n); |
| globpr=1; /* to print the contributions */ | dcwave=ivector(1,n); |
| likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ | |
| printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); | for (i=1; i<=imx; i++){ |
| for (k=1; k<=npar;k++) | dcwave[i]=-1; |
| printf(" %d %8.5f",k,p[k]); | for (j=1; j<=lastpass; j++) |
| printf("\n"); | if (s[j][i]>nlstate) { |
| if(mle>=1){ /* Could be 1 or 2 */ | dcwave[i]=j; |
| mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); | /* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/ |
| break; | |
| } | |
| } | |
| for (i=1; i<=imx; i++) { | |
| if (wav[i]>0){ | |
| ageexmed[i]=agev[mw[1][i]][i]; | |
| j=wav[i];agecens[i]=1.; | |
| if (ageexmed[i]>1 & wav[i]>0) agecens[i]=agev[mw[j][i]][i]; | |
| cens[i]=1; | |
| if (ageexmed[i]<1) cens[i]=-1; | |
| if (agedc[i]< AGESUP & agedc[i]>1 & dcwave[i]>firstpass & dcwave[i]<=lastpass) cens[i]=0 ; | |
| } | |
| else cens[i]=-1; | |
| } | |
| for (i=1;i<=NDIM;i++) { | |
| for (j=1;j<=NDIM;j++) | |
| ximort[i][j]=(i == j ? 1.0 : 0.0); | |
| } | |
| p[1]=0.1; p[2]=0.1; | |
| /*printf("%lf %lf", p[1], p[2]);*/ | |
| printf("Powell\n"); fprintf(ficlog,"Powell\n"); | |
| strcpy(filerespow,"pow-mort"); | |
| strcat(filerespow,fileres); | |
| if((ficrespow=fopen(filerespow,"w"))==NULL) { | |
| printf("Problem with resultfile: %s\n", filerespow); | |
| fprintf(ficlog,"Problem with resultfile: %s\n", filerespow); | |
| } | } |
| fprintf(ficrespow,"# Powell\n# iter -2*LL"); | |
| /* for (i=1;i<=nlstate;i++) | |
| for(j=1;j<=nlstate+ndeath;j++) | |
| if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); | |
| */ | |
| fprintf(ficrespow,"\n"); | |
| powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); | |
| fclose(ficrespow); | |
| /*--------- results files --------------*/ | hesscov(matcov, p, NDIM,delti, 1e-4, gompertz); |
| fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model); | |
| jk=1; | for(i=1; i <=NDIM; i++) |
| fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); | for(j=i+1;j<=NDIM;j++) |
| printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); | matcov[i][j]=matcov[j][i]; |
| fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); | |
| for(i=1,jk=1; i <=nlstate; i++){ | printf("\nCovariance matrix\n "); |
| for(k=1; k <=(nlstate+ndeath); k++){ | for(i=1; i <=NDIM; i++) { |
| if (k != i) | for(j=1;j<=NDIM;j++){ |
| { | printf("%f ",matcov[i][j]); |
| } | |
| printf("\n "); | |
| } | |
| printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp); | |
| for (i=1;i<=NDIM;i++) | |
| printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); | |
| lsurv=vector(1,AGESUP); | |
| lpop=vector(1,AGESUP); | |
| tpop=vector(1,AGESUP); | |
| lsurv[agegomp]=100000; | |
| for (k=agegomp;k<=AGESUP;k++) { | |
| agemortsup=k; | |
| if (p[1]*exp(p[2]*(k-agegomp))>1) break; | |
| } | |
| for (k=agegomp;k<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 / */ | |
| printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); | |
| printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \ | |
| stepm, weightopt,\ | |
| model,imx,p,matcov,agemortsup); | |
| free_vector(lsurv,1,AGESUP); | |
| free_vector(lpop,1,AGESUP); | |
| free_vector(tpop,1,AGESUP); | |
| } /* Endof if mle==-3 */ | |
| else{ /* For mle >=1 */ | |
| likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ | |
| printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); | |
| for (k=1; k<=npar;k++) | |
| printf(" %d %8.5f",k,p[k]); | |
| printf("\n"); | |
| globpr=1; /* to print the contributions */ | |
| likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ | |
| printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); | |
| for (k=1; k<=npar;k++) | |
| printf(" %d %8.5f",k,p[k]); | |
| printf("\n"); | |
| if(mle>=1){ /* Could be 1 or 2 */ | |
| mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); | |
| } | |
| /*--------- results files --------------*/ | |
| fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model); | |
| fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); | |
| printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); | |
| fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); | |
| for(i=1,jk=1; i <=nlstate; i++){ | |
| for(k=1; k <=(nlstate+ndeath); k++){ | |
| if (k != i) { | |
| printf("%d%d ",i,k); | printf("%d%d ",i,k); |
| fprintf(ficlog,"%d%d ",i,k); | fprintf(ficlog,"%d%d ",i,k); |
| fprintf(ficres,"%1d%1d ",i,k); | fprintf(ficres,"%1d%1d ",i,k); |
| Line 4641 Interval (in months) between two waves: | Line 4971 Interval (in months) between two waves: |
| fprintf(ficlog,"\n"); | fprintf(ficlog,"\n"); |
| fprintf(ficres,"\n"); | fprintf(ficres,"\n"); |
| } | } |
| } | |
| } | } |
| } | if(mle!=0){ |
| if(mle!=0){ | /* Computing hessian and covariance matrix */ |
| /* Computing hessian and covariance matrix */ | ftolhess=ftol; /* Usually correct */ |
| ftolhess=ftol; /* Usually correct */ | hesscov(matcov, p, npar, delti, ftolhess, func); |
| hesscov(matcov, p, npar, delti, ftolhess, func); | } |
| } | fprintf(ficres,"# Scales (for hessian or gradient estimation)\n"); |
| fprintf(ficres,"# Scales (for hessian or gradient estimation)\n"); | printf("# Scales (for hessian or gradient estimation)\n"); |
| printf("# Scales (for hessian or gradient estimation)\n"); | fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n"); |
| fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n"); | for(i=1,jk=1; i <=nlstate; i++){ |
| for(i=1,jk=1; i <=nlstate; i++){ | for(j=1; j <=nlstate+ndeath; j++){ |
| for(j=1; j <=nlstate+ndeath; j++){ | if (j!=i) { |
| if (j!=i) { | fprintf(ficres,"%1d%1d",i,j); |
| fprintf(ficres,"%1d%1d",i,j); | printf("%1d%1d",i,j); |
| printf("%1d%1d",i,j); | fprintf(ficlog,"%1d%1d",i,j); |
| fprintf(ficlog,"%1d%1d",i,j); | for(k=1; k<=ncovmodel;k++){ |
| for(k=1; k<=ncovmodel;k++){ | printf(" %.5e",delti[jk]); |
| printf(" %.5e",delti[jk]); | fprintf(ficlog," %.5e",delti[jk]); |
| fprintf(ficlog," %.5e",delti[jk]); | fprintf(ficres," %.5e",delti[jk]); |
| fprintf(ficres," %.5e",delti[jk]); | jk++; |
| jk++; | } |
| printf("\n"); | |
| fprintf(ficlog,"\n"); | |
| fprintf(ficres,"\n"); | |
| } | } |
| printf("\n"); | |
| fprintf(ficlog,"\n"); | |
| fprintf(ficres,"\n"); | |
| } | } |
| } | } |
| } | |
| fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); | |
| fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); | if(mle>=1) |
| if(mle==1) | printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); |
| printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); | fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); |
| fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); | /* # 121 Var(a12)\n\ */ |
| for(i=1,k=1;i<=npar;i++){ | /* # 122 Cov(b12,a12) Var(b12)\n\ */ |
| /* if (k>nlstate) k=1; | /* # 131 Cov(a13,a12) Cov(a13,b12, Var(a13)\n\ */ |
| i1=(i-1)/(ncovmodel*nlstate)+1; | /* # 132 Cov(b13,a12) Cov(b13,b12, Cov(b13,a13) Var(b13)\n\ */ |
| fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]); | /* # 212 Cov(a21,a12) Cov(a21,b12, Cov(a21,a13) Cov(a21,b13) Var(a21)\n\ */ |
| printf("%s%d%d",alph[k],i1,tab[i]); | /* # 212 Cov(b21,a12) Cov(b21,b12, Cov(b21,a13) Cov(b21,b13) Cov(b21,a21) Var(b21)\n\ */ |
| /* # 232 Cov(a23,a12) Cov(a23,b12, Cov(a23,a13) Cov(a23,b13) Cov(a23,a21) Cov(a23,b21) Var(a23)\n\ */ | |
| /* # 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n" */ | |
| /* Just to have a covariance matrix which will be more understandable | |
| even is we still don't want to manage dictionary of variables | |
| */ | */ |
| fprintf(ficres,"%3d",i); | for(itimes=1;itimes<=2;itimes++){ |
| if(mle==1) | jj=0; |
| printf("%3d",i); | for(i=1; i <=nlstate; i++){ |
| fprintf(ficlog,"%3d",i); | for(j=1; j <=nlstate+ndeath; j++){ |
| for(j=1; j<=i;j++){ | if(j==i) continue; |
| fprintf(ficres," %.5e",matcov[i][j]); | for(k=1; k<=ncovmodel;k++){ |
| if(mle==1) | jj++; |
| printf(" %.5e",matcov[i][j]); | ca[0]= k+'a'-1;ca[1]='\0'; |
| fprintf(ficlog," %.5e",matcov[i][j]); | if(itimes==1){ |
| if(mle>=1) | |
| printf("#%1d%1d%d",i,j,k); | |
| fprintf(ficlog,"#%1d%1d%d",i,j,k); | |
| fprintf(ficres,"#%1d%1d%d",i,j,k); | |
| }else{ | |
| if(mle>=1) | |
| printf("%1d%1d%d",i,j,k); | |
| fprintf(ficlog,"%1d%1d%d",i,j,k); | |
| fprintf(ficres,"%1d%1d%d",i,j,k); | |
| } | |
| ll=0; | |
| for(li=1;li <=nlstate; li++){ | |
| for(lj=1;lj <=nlstate+ndeath; lj++){ | |
| if(lj==li) continue; | |
| for(lk=1;lk<=ncovmodel;lk++){ | |
| ll++; | |
| if(ll<=jj){ | |
| cb[0]= lk +'a'-1;cb[1]='\0'; | |
| if(ll<jj){ | |
| if(itimes==1){ | |
| if(mle>=1) | |
| printf(" Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj); | |
| fprintf(ficlog," Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj); | |
| fprintf(ficres," Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj); | |
| }else{ | |
| if(mle>=1) | |
| printf(" %.5e",matcov[jj][ll]); | |
| fprintf(ficlog," %.5e",matcov[jj][ll]); | |
| fprintf(ficres," %.5e",matcov[jj][ll]); | |
| } | |
| }else{ | |
| if(itimes==1){ | |
| if(mle>=1) | |
| printf(" Var(%s%1d%1d)",ca,i,j); | |
| fprintf(ficlog," Var(%s%1d%1d)",ca,i,j); | |
| fprintf(ficres," Var(%s%1d%1d)",ca,i,j); | |
| }else{ | |
| if(mle>=1) | |
| printf(" %.5e",matcov[jj][ll]); | |
| fprintf(ficlog," %.5e",matcov[jj][ll]); | |
| fprintf(ficres," %.5e",matcov[jj][ll]); | |
| } | |
| } | |
| } | |
| } /* end lk */ | |
| } /* end lj */ | |
| } /* end li */ | |
| if(mle>=1) | |
| printf("\n"); | |
| fprintf(ficlog,"\n"); | |
| fprintf(ficres,"\n"); | |
| numlinepar++; | |
| } /* end k*/ | |
| } /*end j */ | |
| } /* end i */ | |
| } /* end itimes */ | |
| fflush(ficlog); | |
| fflush(ficres); | |
| while((c=getc(ficpar))=='#' && c!= EOF){ | |
| ungetc(c,ficpar); | |
| fgets(line, MAXLINE, ficpar); | |
| puts(line); | |
| fputs(line,ficparo); | |
| } | } |
| fprintf(ficres,"\n"); | |
| if(mle==1) | |
| printf("\n"); | |
| fprintf(ficlog,"\n"); | |
| k++; | |
| } | |
| while((c=getc(ficpar))=='#' && c!= EOF){ | |
| ungetc(c,ficpar); | |
| fgets(line, MAXLINE, ficpar); | |
| puts(line); | |
| fputs(line,ficparo); | |
| } | |
| ungetc(c,ficpar); | |
| estepm=0; | |
| fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); | |
| if (estepm==0 || estepm < stepm) estepm=stepm; | |
| if (fage <= 2) { | |
| bage = ageminpar; | |
| fage = agemaxpar; | |
| } | |
| fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); | |
| fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); | |
| fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); | |
| while((c=getc(ficpar))=='#' && c!= EOF){ | |
| ungetc(c,ficpar); | ungetc(c,ficpar); |
| fgets(line, MAXLINE, ficpar); | |
| puts(line); | estepm=0; |
| fputs(line,ficparo); | fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); |
| } | if (estepm==0 || estepm < stepm) estepm=stepm; |
| ungetc(c,ficpar); | if (fage <= 2) { |
| bage = ageminpar; | |
| fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav); | fage = agemaxpar; |
| fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); | } |
| fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); | |
| printf("begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); | fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); |
| fprintf(ficlog,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); | fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); |
| fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); | |
| while((c=getc(ficpar))=='#' && c!= EOF){ | |
| while((c=getc(ficpar))=='#' && c!= EOF){ | |
| ungetc(c,ficpar); | |
| fgets(line, MAXLINE, ficpar); | |
| puts(line); | |
| fputs(line,ficparo); | |
| } | |
| ungetc(c,ficpar); | ungetc(c,ficpar); |
| fgets(line, MAXLINE, ficpar); | |
| puts(line); | fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav); |
| fputs(line,ficparo); | fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); |
| } | fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); |
| ungetc(c,ficpar); | printf("begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); |
| fprintf(ficlog,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); | |
| dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.; | while((c=getc(ficpar))=='#' && c!= EOF){ |
| dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.; | ungetc(c,ficpar); |
| fgets(line, MAXLINE, ficpar); | |
| fscanf(ficpar,"pop_based=%d\n",&popbased); | puts(line); |
| fprintf(ficparo,"pop_based=%d\n",popbased); | fputs(line,ficparo); |
| fprintf(ficres,"pop_based=%d\n",popbased); | } |
| while((c=getc(ficpar))=='#' && c!= EOF){ | |
| ungetc(c,ficpar); | ungetc(c,ficpar); |
| fgets(line, MAXLINE, ficpar); | |
| puts(line); | |
| fputs(line,ficparo); | dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.; |
| } | dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.; |
| ungetc(c,ficpar); | |
| fscanf(ficpar,"pop_based=%d\n",&popbased); | |
| fscanf(ficpar,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj); | fprintf(ficparo,"pop_based=%d\n",popbased); |
| fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); | fprintf(ficres,"pop_based=%d\n",popbased); |
| printf("prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); | |
| fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); | while((c=getc(ficpar))=='#' && c!= EOF){ |
| fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); | ungetc(c,ficpar); |
| /* day and month of proj2 are not used but only year anproj2.*/ | fgets(line, MAXLINE, ficpar); |
| puts(line); | |
| while((c=getc(ficpar))=='#' && c!= EOF){ | fputs(line,ficparo); |
| } | |
| ungetc(c,ficpar); | ungetc(c,ficpar); |
| fgets(line, MAXLINE, ficpar); | |
| puts(line); | fscanf(ficpar,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj); |
| fputs(line,ficparo); | fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); |
| } | printf("prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); |
| ungetc(c,ficpar); | fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); |
| fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); | |
| fscanf(ficpar,"popforecast=%d popfile=%s popfiledate=%lf/%lf/%lf last-popfiledate=%lf/%lf/%lf\n",&popforecast,popfile,&jpyram,&mpyram,&anpyram,&jpyram1,&mpyram1,&anpyram1); | /* day and month of proj2 are not used but only year anproj2.*/ |
| fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); | |
| fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); | |
| /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);*/ | /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);*/ |
| /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ | /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
| replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */ | replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */ |
| printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); | printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); |
| printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\ | printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\ |
| model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\ | model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\ |
| jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); | jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); |
| /*------------ free_vector -------------*/ | /*------------ free_vector -------------*/ |
| /* chdir(path); */ | /* chdir(path); */ |
| free_ivector(wav,1,imx); | free_ivector(wav,1,imx); |
| free_imatrix(dh,1,lastpass-firstpass+1,1,imx); | free_imatrix(dh,1,lastpass-firstpass+1,1,imx); |
| free_imatrix(bh,1,lastpass-firstpass+1,1,imx); | free_imatrix(bh,1,lastpass-firstpass+1,1,imx); |
| free_imatrix(mw,1,lastpass-firstpass+1,1,imx); | free_imatrix(mw,1,lastpass-firstpass+1,1,imx); |
| free_lvector(num,1,n); | free_lvector(num,1,n); |
| free_vector(agedc,1,n); | free_vector(agedc,1,n); |
| /*free_matrix(covar,0,NCOVMAX,1,n);*/ | /*free_matrix(covar,0,NCOVMAX,1,n);*/ |
| /*free_matrix(covar,1,NCOVMAX,1,n);*/ | /*free_matrix(covar,1,NCOVMAX,1,n);*/ |
| fclose(ficparo); | fclose(ficparo); |
| fclose(ficres); | fclose(ficres); |
| /*--------------- Prevalence limit (stable prevalence) --------------*/ | /*--------------- Prevalence limit (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 stable prevalence resultfile: %s\n", filerespl);goto end; |
| fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end; | fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end; |
| } | } |
| printf("Computing stable prevalence: result on file '%s' \n", filerespl); | printf("Computing stable prevalence: result on file '%s' \n", filerespl); |
| fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl); | fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl); |
| fprintf(ficrespl,"#Stable prevalence \n"); | fprintf(ficrespl, "#Local time at start: %s", strstart); |
| fprintf(ficrespl,"#Age "); | fprintf(ficrespl,"#Stable prevalence \n"); |
| for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); | fprintf(ficrespl,"#Age "); |
| fprintf(ficrespl,"\n"); | for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); |
| fprintf(ficrespl,"\n"); | |
| prlim=matrix(1,nlstate,1,nlstate); | |
| prlim=matrix(1,nlstate,1,nlstate); | |
| agebase=ageminpar; | |
| agelim=agemaxpar; | agebase=ageminpar; |
| ftolpl=1.e-10; | agelim=agemaxpar; |
| i1=cptcoveff; | ftolpl=1.e-10; |
| if (cptcovn < 1){i1=1;} | i1=cptcoveff; |
| if (cptcovn < 1){i1=1;} | |
| for(cptcov=1,k=0;cptcov<=i1;cptcov++){ | |
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ | for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| k=k+1; | for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
| /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/ | k=k+1; |
| fprintf(ficrespl,"\n#******"); | /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/ |
| printf("\n#******"); | fprintf(ficrespl,"\n#******"); |
| fprintf(ficlog,"\n#******"); | printf("\n#******"); |
| for(j=1;j<=cptcoveff;j++) { | fprintf(ficlog,"\n#******"); |
| fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); | for(j=1;j<=cptcoveff;j++) { |
| printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); | fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
| fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); | printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
| } | fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
| fprintf(ficrespl,"******\n"); | } |
| printf("******\n"); | fprintf(ficrespl,"******\n"); |
| fprintf(ficlog,"******\n"); | printf("******\n"); |
| fprintf(ficlog,"******\n"); | |
| for (age=agebase; age<=agelim; age++){ | for (age=agebase; age<=agelim; age++){ |
| prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); | prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); |
| fprintf(ficrespl,"%.0f ",age ); | fprintf(ficrespl,"%.0f ",age ); |
| for(j=1;j<=cptcoveff;j++) | for(j=1;j<=cptcoveff;j++) |
| fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); | fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
| for(i=1; i<=nlstate;i++) | for(i=1; i<=nlstate;i++) |
| fprintf(ficrespl," %.5f", prlim[i][i]); | fprintf(ficrespl," %.5f", prlim[i][i]); |
| fprintf(ficrespl,"\n"); | fprintf(ficrespl,"\n"); |
| } | |
| } | } |
| } | } |
| } | fclose(ficrespl); |
| fclose(ficrespl); | |
| /*------------- h Pij x at various ages ------------*/ | /*------------- h Pij x at various ages ------------*/ |
| strcpy(filerespij,"pij"); strcat(filerespij,fileres); | strcpy(filerespij,"pij"); strcat(filerespij,fileres); |
| if((ficrespij=fopen(filerespij,"w"))==NULL) { | if((ficrespij=fopen(filerespij,"w"))==NULL) { |
| printf("Problem with Pij resultfile: %s\n", filerespij);goto end; | printf("Problem with Pij resultfile: %s\n", filerespij);goto end; |
| fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij);goto end; | fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij);goto end; |
| } | } |
| printf("Computing pij: result on file '%s' \n", filerespij); | printf("Computing pij: result on file '%s' \n", filerespij); |
| fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij); | fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij); |
| stepsize=(int) (stepm+YEARM-1)/YEARM; | stepsize=(int) (stepm+YEARM-1)/YEARM; |
| /*if (stepm<=24) stepsize=2;*/ | /*if (stepm<=24) stepsize=2;*/ |
| agelim=AGESUP; | |
| hstepm=stepsize*YEARM; /* Every year of age */ | |
| hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ | |
| /* hstepm=1; aff par mois*/ | agelim=AGESUP; |
| hstepm=stepsize*YEARM; /* Every year of age */ | |
| hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ | |
| fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x "); | /* hstepm=1; aff par mois*/ |
| for(cptcov=1,k=0;cptcov<=i1;cptcov++){ | fprintf(ficrespij, "#Local time at start: %s", strstart); |
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ | fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x "); |
| k=k+1; | for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| fprintf(ficrespij,"\n#****** "); | for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
| for(j=1;j<=cptcoveff;j++) | k=k+1; |
| fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); | fprintf(ficrespij,"\n#****** "); |
| fprintf(ficrespij,"******\n"); | for(j=1;j<=cptcoveff;j++) |
| fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); | |
| fprintf(ficrespij,"******\n"); | |
| for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */ | for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */ |
| nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ | nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ |
| nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ | nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ |
| /* nhstepm=nhstepm*YEARM; aff par mois*/ | /* nhstepm=nhstepm*YEARM; aff par mois*/ |
| p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); | p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
| oldm=oldms;savm=savms; | oldm=oldms;savm=savms; |
| hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); | hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); |
| fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j="); | fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j="); |
| for(i=1; i<=nlstate;i++) | |
| for(j=1; j<=nlstate+ndeath;j++) | |
| fprintf(ficrespij," %1d-%1d",i,j); | |
| fprintf(ficrespij,"\n"); | |
| for (h=0; h<=nhstepm; h++){ | |
| fprintf(ficrespij,"%d %3.f %3.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm ); | |
| 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++) |
| fprintf(ficrespij," %.5f", p3mat[i][j][h]); | fprintf(ficrespij," %1d-%1d",i,j); |
| fprintf(ficrespij,"\n"); | |
| for (h=0; h<=nhstepm; h++){ | |
| fprintf(ficrespij,"%d %3.f %3.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm ); | |
| for(i=1; i<=nlstate;i++) | |
| for(j=1; j<=nlstate+ndeath;j++) | |
| fprintf(ficrespij," %.5f", p3mat[i][j][h]); | |
| fprintf(ficrespij,"\n"); | |
| } | |
| free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); | |
| fprintf(ficrespij,"\n"); | fprintf(ficrespij,"\n"); |
| } | } |
| free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); | |
| fprintf(ficrespij,"\n"); | |
| } | } |
| } | } |
| } | |
| varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax); | |
| fclose(ficrespij); | varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart); |
| probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); | fclose(ficrespij); |
| /*---------- Forecasting ------------------*/ | probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| /*if((stepm == 1) && (strcmp(model,".")==0)){*/ | for(i=1;i<=AGESUP;i++) |
| if(prevfcast==1){ | for(j=1;j<=NCOVMAX;j++) |
| /* if(stepm ==1){*/ | for(k=1;k<=NCOVMAX;k++) |
| probs[i][j][k]=0.; | |
| /*---------- Forecasting ------------------*/ | |
| /*if((stepm == 1) && (strcmp(model,".")==0)){*/ | |
| if(prevfcast==1){ | |
| /* if(stepm ==1){*/ | |
| prevforecast(fileres, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff); | prevforecast(fileres, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff); |
| /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/ | /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/ |
| /* } */ | /* } */ |
| /* else{ */ | /* else{ */ |
| /* erreur=108; */ | /* erreur=108; */ |
| /* printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */ | /* printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */ |
| /* fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */ | /* fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */ |
| /* } */ | /* } */ |
| } | } |
| /*---------- Health expectancies and variances ------------*/ | /*---------- Health expectancies and variances ------------*/ |
| strcpy(filerest,"t"); | strcpy(filerest,"t"); |
| strcat(filerest,fileres); | strcat(filerest,fileres); |
| if((ficrest=fopen(filerest,"w"))==NULL) { | if((ficrest=fopen(filerest,"w"))==NULL) { |
| 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 LEs with variances: file '%s' \n", filerest); |
| fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest); | fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest); |
| strcpy(filerese,"e"); | strcpy(filerese,"e"); |
| strcat(filerese,fileres); | strcat(filerese,fileres); |
| if((ficreseij=fopen(filerese,"w"))==NULL) { | if((ficreseij=fopen(filerese,"w"))==NULL) { |
| printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0); | printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0); |
| fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0); | fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0); |
| } | } |
| printf("Computing Health Expectancies: result on file '%s' \n", filerese); | 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(fileresv,"v"); | strcpy(fileresv,"v"); |
| strcat(fileresv,fileres); | strcat(fileresv,fileres); |
| if((ficresvij=fopen(fileresv,"w"))==NULL) { | if((ficresvij=fopen(fileresv,"w"))==NULL) { |
| printf("Problem with variance resultfile: %s\n", fileresv);exit(0); | printf("Problem with variance resultfile: %s\n", fileresv);exit(0); |
| fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0); | fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0); |
| } | } |
| printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); | printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); |
| fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); | fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); |
| /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */ | /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */ |
| prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); | prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
| /* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\ | /* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\ |
| ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass); | ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass); |
| */ | */ |
| if (mobilav!=0) { | if (mobilav!=0) { |
| mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); | mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ | if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ |
| fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); | fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); |
| printf(" Error in movingaverage mobilav=%d\n",mobilav); | printf(" Error in movingaverage mobilav=%d\n",mobilav); |
| } | |
| } | } |
| } | |
| 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++){ |
| k=k+1; | k=k+1; |
| fprintf(ficrest,"\n#****** "); | fprintf(ficrest,"\n#****** "); |
| for(j=1;j<=cptcoveff;j++) | for(j=1;j<=cptcoveff;j++) |
| fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); | fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
| fprintf(ficrest,"******\n"); | fprintf(ficrest,"******\n"); |
| fprintf(ficreseij,"\n#****** "); | fprintf(ficreseij,"\n#****** "); |
| for(j=1;j<=cptcoveff;j++) | 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(ficreseij,"******\n"); | fprintf(ficreseij,"******\n"); |
| fprintf(ficresvij,"\n#****** "); | fprintf(ficresvij,"\n#****** "); |
| for(j=1;j<=cptcoveff;j++) | for(j=1;j<=cptcoveff;j++) |
| fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); | fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
| fprintf(ficresvij,"******\n"); | fprintf(ficresvij,"******\n"); |
| eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); | |
| oldm=oldms;savm=savms; | |
| evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov); | |
| vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); | |
| oldm=oldms;savm=savms; | |
| varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav); | |
| if(popbased==1){ | |
| varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav); | |
| } | |
| eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); | |
| oldm=oldms;savm=savms; | |
| evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart); | |
| fprintf(ficrest,"#Total LEs with variances: e.. (std) "); | vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); |
| for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); | oldm=oldms;savm=savms; |
| fprintf(ficrest,"\n"); | 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){ | |
| epj=vector(1,nlstate+1); | varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav, strstart); |
| for(age=bage; age <=fage ;age++){ | |
| prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); | |
| if (popbased==1) { | |
| if(mobilav ==0){ | |
| for(i=1; i<=nlstate;i++) | |
| prlim[i][i]=probs[(int)age][i][k]; | |
| }else{ /* mobilav */ | |
| for(i=1; i<=nlstate;i++) | |
| prlim[i][i]=mobaverage[(int)age][i][k]; | |
| } | |
| } | } |
| fprintf(ficrest, "#Local time at start: %s", strstart); | |
| fprintf(ficrest,"#Total LEs with variances: e.. (std) "); | |
| for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); | |
| fprintf(ficrest,"\n"); | |
| epj=vector(1,nlstate+1); | |
| for(age=bage; age <=fage ;age++){ | |
| prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); | |
| if (popbased==1) { | |
| if(mobilav ==0){ | |
| for(i=1; i<=nlstate;i++) | |
| prlim[i][i]=probs[(int)age][i][k]; | |
| }else{ /* mobilav */ | |
| for(i=1; i<=nlstate;i++) | |
| prlim[i][i]=mobaverage[(int)age][i][k]; | |
| } | |
| } | |
| fprintf(ficrest," %4.0f",age); | fprintf(ficrest," %4.0f",age); |
| for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ | for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ |
| for(i=1, epj[j]=0.;i <=nlstate;i++) { | for(i=1, epj[j]=0.;i <=nlstate;i++) { |
| epj[j] += prlim[i][i]*eij[i][j][(int)age]; | epj[j] += prlim[i][i]*eij[i][j][(int)age]; |
| /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ | /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ |
| } | |
| epj[nlstate+1] +=epj[j]; | |
| } | } |
| epj[nlstate+1] +=epj[j]; | |
| } | |
| for(i=1, vepp=0.;i <=nlstate;i++) | for(i=1, vepp=0.;i <=nlstate;i++) |
| for(j=1;j <=nlstate;j++) | for(j=1;j <=nlstate;j++) |
| vepp += vareij[i][j][(int)age]; | vepp += vareij[i][j][(int)age]; |
| fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); | fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); |
| for(j=1;j <=nlstate;j++){ | for(j=1;j <=nlstate;j++){ |
| fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); | fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); |
| } | |
| fprintf(ficrest,"\n"); | |
| } | } |
| fprintf(ficrest,"\n"); | free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
| } | free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
| free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); | free_vector(epj,1,nlstate+1); |
| free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); | } |
| free_vector(epj,1,nlstate+1); | } |
| } | free_vector(weight,1,n); |
| } | free_imatrix(Tvard,1,15,1,2); |
| free_vector(weight,1,n); | free_imatrix(s,1,maxwav+1,1,n); |
| free_imatrix(Tvard,1,15,1,2); | free_matrix(anint,1,maxwav,1,n); |
| free_imatrix(s,1,maxwav+1,1,n); | free_matrix(mint,1,maxwav,1,n); |
| free_matrix(anint,1,maxwav,1,n); | free_ivector(cod,1,n); |
| free_matrix(mint,1,maxwav,1,n); | free_ivector(tab,1,NCOVMAX); |
| free_ivector(cod,1,n); | fclose(ficreseij); |
| free_ivector(tab,1,NCOVMAX); | fclose(ficresvij); |
| fclose(ficreseij); | fclose(ficrest); |
| fclose(ficresvij); | fclose(ficpar); |
| fclose(ficrest); | |
| fclose(ficpar); | /*------- Variance of stable prevalence------*/ |
| /*------- Variance of stable prevalence------*/ | strcpy(fileresvpl,"vpl"); |
| strcat(fileresvpl,fileres); | |
| strcpy(fileresvpl,"vpl"); | if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { |
| strcat(fileresvpl,fileres); | printf("Problem with variance of stable prevalence resultfile: %s\n", fileresvpl); |
| if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { | exit(0); |
| printf("Problem with variance of stable prevalence resultfile: %s\n", fileresvpl); | } |
| exit(0); | printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl); |
| } | |
| printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl); | for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ | |
| for(cptcov=1,k=0;cptcov<=i1;cptcov++){ | k=k+1; |
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ | fprintf(ficresvpl,"\n#****** "); |
| k=k+1; | for(j=1;j<=cptcoveff;j++) |
| fprintf(ficresvpl,"\n#****** "); | fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
| for(j=1;j<=cptcoveff;j++) | fprintf(ficresvpl,"******\n"); |
| fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); | |
| fprintf(ficresvpl,"******\n"); | |
| 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); |
| } | |
| } | } |
| } | |
| fclose(ficresvpl); | fclose(ficresvpl); |
| /*---------- End : free ----------------*/ | /*---------- End : free ----------------*/ |
| free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); | if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); | free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); | |
| free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); | } /* mle==-3 arrives here for freeing */ |
| free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); | |
| free_matrix(covar,0,NCOVMAX,1,n); | free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(matcov,1,npar,1,npar); | free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| /*free_vector(delti,1,npar);*/ | free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); | |
| free_matrix(agev,1,maxwav,1,imx); | free_matrix(covar,0,NCOVMAX,1,n); |
| free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); | free_matrix(matcov,1,npar,1,npar); |
| if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); | /*free_vector(delti,1,npar);*/ |
| free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); | free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); |
| free_matrix(agev,1,maxwav,1,imx); | |
| free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); | |
| free_ivector(ncodemax,1,8); | |
| free_ivector(Tvar,1,15); | |
| free_ivector(Tprod,1,15); | |
| free_ivector(Tvaraff,1,15); | |
| free_ivector(Tage,1,15); | |
| free_ivector(Tcode,1,100); | |
| free_ivector(ncodemax,1,8); | |
| free_ivector(Tvar,1,15); | |
| free_ivector(Tprod,1,15); | |
| free_ivector(Tvaraff,1,15); | |
| free_ivector(Tage,1,15); | |
| free_ivector(Tcode,1,100); | |
| fflush(fichtm); | fflush(fichtm); |
| fflush(ficgp); | fflush(ficgp); |
| if((nberr >0) || (nbwarn>0)){ | if((nberr >0) || (nbwarn>0)){ |
| printf("End of Imach with %d errors and/or warnings %d\n",nberr,nbwarn); | printf("End of Imach with %d errors and/or %d warnings\n",nberr,nbwarn); |
| fprintf(ficlog,"End of Imach with %d errors and/or warnings %d\n",nberr,nbwarn); | fprintf(ficlog,"End of Imach with %d errors and/or warnings %d\n",nberr,nbwarn); |
| }else{ | }else{ |
| printf("End of Imach\n"); | printf("End of Imach\n"); |
| Line 5119 ageminpar, agemax, s[lastpass][imx], age | Line 5506 ageminpar, agemax, s[lastpass][imx], age |
| 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("Localtime at start %s\nLocaltime at end %s",strstart, strtend); | printf("Local time at start %s\nLocal time at end %s",strstart, strtend); |
| fprintf(ficlog,"Localtime 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)); |
| printf("Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec); | printf("Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec); |
| Line 5136 ageminpar, agemax, s[lastpass][imx], age | Line 5523 ageminpar, agemax, s[lastpass][imx], age |
| /*------ End -----------*/ | /*------ End -----------*/ |
| chdir(path); | chdir(path); |
| strcpy(plotcmd,GNUPLOTPROGRAM); | strcpy(plotcmd,"\""); |
| strcat(plotcmd,pathimach); | |
| strcat(plotcmd,GNUPLOTPROGRAM); | |
| strcat(plotcmd,"\""); | |
| strcat(plotcmd," "); | strcat(plotcmd," "); |
| strcat(plotcmd,optionfilegnuplot); | strcat(plotcmd,optionfilegnuplot); |
| printf("Starting graphs with: %s",plotcmd);fflush(stdout); | printf("Starting graphs with: %s",plotcmd);fflush(stdout); |
| Line 5149 ageminpar, agemax, s[lastpass][imx], age | Line 5539 ageminpar, agemax, s[lastpass][imx], age |
| printf("\nType e to edit output files, g to graph again and q for exiting: "); | printf("\nType e to edit output files, g to graph again and q for exiting: "); |
| scanf("%s",z); | scanf("%s",z); |
| /* if (z[0] == 'c') system("./imach"); */ | /* if (z[0] == 'c') system("./imach"); */ |
| if (z[0] == 'e') system(optionfilehtm); | if (z[0] == 'e') { |
| printf("Starting browser with: %s",optionfilehtm);fflush(stdout); | |
| system(optionfilehtm); | |
| } | |
| else if (z[0] == 'g') system(plotcmd); | else if (z[0] == 'g') system(plotcmd); |
| else if (z[0] == 'q') exit(0); | else if (z[0] == 'q') exit(0); |
| } | } |