--- imach/src/imach.c 2004/05/16 15:05:56 1.98 +++ imach/src/imach.c 2006/01/30 09:55:26 1.112 @@ -1,6 +1,57 @@ -/* $Id: imach.c,v 1.98 2004/05/16 15:05:56 brouard Exp $ +/* $Id: imach.c,v 1.112 2006/01/30 09:55:26 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.112 2006/01/30 09:55:26 brouard + (Module): Back to gnuplot.exe instead of wgnuplot.exe + + Revision 1.111 2006/01/25 20:38:18 brouard + (Module): Lots of cleaning and bugs added (Gompertz) + (Module): Comments can be added in data file. Missing date values + can be a simple dot '.'. + + Revision 1.110 2006/01/25 00:51:50 brouard + (Module): Lots of cleaning and bugs added (Gompertz) + + Revision 1.109 2006/01/24 19:37:15 brouard + (Module): Comments (lines starting with a #) are allowed in data. + + Revision 1.108 2006/01/19 18:05:42 lievre + Gnuplot problem appeared... + To be fixed + + Revision 1.107 2006/01/19 16:20:37 brouard + Test existence of gnuplot in imach path + + Revision 1.106 2006/01/19 13:24:36 brouard + Some cleaning and links added in html output + + Revision 1.105 2006/01/05 20:23:19 lievre + *** empty log message *** + + Revision 1.104 2005/09/30 16:11:43 lievre + (Module): sump fixed, loop imx fixed, and simplifications. + (Module): If the status is missing at the last wave but we know + that the person is alive, then we can code his/her status as -2 + (instead of missing=-1 in earlier versions) and his/her + contributions to the likelihood is 1 - Prob of dying from last + health status (= 1-p13= p11+p12 in the easiest case of somebody in + the healthy state at last known wave). Version is 0.98 + + Revision 1.103 2005/09/30 15:54:49 lievre + (Module): sump fixed, loop imx fixed, and simplifications. + + Revision 1.102 2004/09/15 17:31:30 brouard + Add the possibility to read data file including tab characters. + + Revision 1.101 2004/09/15 10:38:38 brouard + Fix on curr_time + + Revision 1.100 2004/07/12 18:29:06 brouard + Add version for Mac OS X. Just define UNIX in Makefile + + Revision 1.99 2004/06/05 08:57:40 brouard + *** empty log message *** + Revision 1.98 2004/05/16 15:05:56 brouard New version 0.97 . First attempt to estimate force of mortality directly from the data i.e. without the need of knowing the health @@ -194,9 +245,16 @@ #include #include #include +#include #include -#include +#include +#include +#include +#include +extern int errno; + +/* #include */ #include #include "timeval.h" @@ -204,11 +262,11 @@ /* #define _(String) gettext (String) */ #define MAXLINE 256 + #define GNUPLOTPROGRAM "gnuplot" /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ #define FILENAMELENGTH 132 -/*#define DEBUG*/ -/*#define windows*/ + #define GLOCK_ERROR_NOPATH -1 /* empty path */ #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ @@ -224,19 +282,21 @@ #define AGESUP 130 #define AGEBASE 40 #define AGEGOMP 10. /* Minimal age for Gompertz adjustment */ -#ifdef unix +#ifdef UNIX #define DIRSEPARATOR '/' +#define CHARSEPARATOR "/" #define ODIRSEPARATOR '\\' #else #define DIRSEPARATOR '\\' +#define CHARSEPARATOR "\\" #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.98 2004/05/16 15:05:56 brouard Exp $ */ +/* $Id: imach.c,v 1.112 2006/01/30 09:55:26 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.70, May 2004, INED-EUROREVES "; -char fullversion[]="$Revision: 1.98 $ $Date: 2004/05/16 15:05:56 $"; +char version[]="Imach version 0.98b, January 2006, INED-EUROREVES "; +char fullversion[]="$Revision: 1.112 $ $Date: 2006/01/30 09:55:26 $"; int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ int nvar; int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; @@ -249,6 +309,7 @@ int popbased=0; int *wav; /* Number of waves for this individuual 0 is possible */ int maxwav; /* Maxim number of waves */ int jmin, jmax; /* min, max spacing between 2 waves */ +int ijmin, ijmax; /* Individuals having jmin and jmax */ int gipmx, gsw; /* Global variables on the number of contributions to the likelihood and the sum of weights (done by funcone)*/ int mle, weightopt; @@ -279,7 +340,7 @@ FILE *ficresvpl; char fileresvpl[FILENAMELENGTH]; char title[MAXLINE]; char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH]; -char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH]; +char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH]; char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH]; char command[FILENAMELENGTH]; int outcmd=0; @@ -301,6 +362,9 @@ long time_value; extern long time(); char strcurr[80], strfor[80]; +char *endptr; +long lval; + #define NR_END 1 #define FREE_ARG char* #define FTOL 1.0e-10 @@ -349,6 +413,7 @@ double *weight; int **s; /* Status */ double *agedc, **covar, idx; int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff; +double *lsurv, *lpop, *tpop; double ftol=FTOL; /* Tolerance for computing Max Likelihood */ double ftolhess; /* Tolerance for computing hessian */ @@ -356,13 +421,17 @@ double ftolhess; /* Tolerance for comput /**************** split *************************/ static int split( char *path, char *dirc, char *name, char *ext, char *finame ) { + /* From a file name with (full) path (either Unix or Windows) we extract the directory (dirc) + the name of the file (name), its extension only (ext) and its first part of the name (finame) + */ char *ss; /* pointer */ int l1, l2; /* length counters */ l1 = strlen(path ); /* length of path */ if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH ); ss= strrchr( path, DIRSEPARATOR ); /* find last / */ - if ( ss == NULL ) { /* no directory, so use current */ + if ( ss == NULL ) { /* no directory, so determine current directory */ + strcpy( name, path ); /* we got the fullname name because no directory */ /*if(strrchr(path, ODIRSEPARATOR )==NULL) printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/ /* get current working directory */ @@ -370,7 +439,8 @@ static int split( char *path, char *dirc if ( getcwd( dirc, FILENAME_MAX ) == NULL ) { return( GLOCK_ERROR_GETCWD ); } - strcpy( name, path ); /* we've got it */ + /* got dirc from getcwd*/ + printf(" DIRC = %s \n",dirc); } else { /* strip direcotry from path */ ss++; /* after this, the filename */ l2 = strlen( ss ); /* length of filename */ @@ -378,21 +448,25 @@ static int split( char *path, char *dirc strcpy( name, ss ); /* save file name */ strncpy( dirc, path, l1 - l2 ); /* now the directory */ dirc[l1-l2] = 0; /* add zero */ + printf(" DIRC2 = %s \n",dirc); } + /* We add a separator at the end of dirc if not exists */ l1 = strlen( dirc ); /* length of directory */ - /*#ifdef windows - if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; } -#else - if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; } -#endif - */ + if( dirc[l1-1] != DIRSEPARATOR ){ + dirc[l1] = DIRSEPARATOR; + dirc[l1+1] = 0; + printf(" DIRC3 = %s \n",dirc); + } ss = strrchr( name, '.' ); /* find last / */ - ss++; - strcpy(ext,ss); /* save extension */ - l1= strlen( name); - l2= strlen(ss)+1; - strncpy( finame, name, l1-l2); - finame[l1-l2]= 0; + if (ss >0){ + ss++; + strcpy(ext,ss); /* save extension */ + l1= strlen( name); + l2= strlen(ss)+1; + strncpy( finame, name, l1-l2); + finame[l1-l2]= 0; + } + return( 0 ); /* we're done */ } @@ -425,8 +499,8 @@ int nbocc(char *s, char occ) void cutv(char *u,char *v, char*t, char occ) { - /* cuts string t into u and v where u is ended by char occ excluding it - and v is after occ excluding it too : ex cutv(u,v,"abcdef2ghi2j",2) + /* cuts string t into u and v where u ends before first occurence of char 'occ' + and v starts after first occurence of char 'occ' : ex cutv(u,v,"abcdef2ghi2j",'2') gives u="abcedf" and v="ghi2j" */ int i,lg,j,p=0; i=0; @@ -868,11 +942,11 @@ void powell(double p[], double **xi, int fprintf(ficrespow,"\n");fflush(ficrespow); if(*iter <=3){ tm = *localtime(&curr_time.tv_sec); - strcpy(strcurr,asctime(&tmf)); + strcpy(strcurr,asctime(&tm)); /* asctime_r(&tm,strcurr); */ - forecast_time=curr_time; + forecast_time=curr_time; itmp = strlen(strcurr); - if(strcurr[itmp-1]=='\n') + if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */ strcurr[itmp-1]='\0'; printf("\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec); fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec); @@ -884,8 +958,8 @@ void powell(double p[], double **xi, int itmp = strlen(strfor); if(strfor[itmp-1]=='\n') strfor[itmp-1]='\0'; - printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s or\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); - fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s or\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); + printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); + fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); } } for (i=1;i<=n;i++) { @@ -1050,57 +1124,57 @@ double **pmij(double **ps, double *cov, int i,j,j1, nc, ii, jj; for(i=1; i<= nlstate; i++){ - for(j=1; ji s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/ + for(j=1; ji s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2); */ + } + ps[i][j]=s2; } - ps[i][j]=s2; } - } /*ps[3][2]=1;*/ - - for(i=1; i<= nlstate; i++){ - s1=0; - for(j=1; j 1 the results are less biased than in previous versions. @@ -1239,15 +1313,16 @@ double func( double *x) s1=s[mw[mi][i]][i]; s2=s[mw[mi+1][i]][i]; bbh=(double)bh[mi][i]/(double)stepm; - /* bias is positive if real duration + /* bias bh is positive if real duration * is higher than the multiple of stepm and negative otherwise. */ /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ if( s2 > nlstate){ - /* i.e. if s2 is a death state and if the date of death is known then the contribution - to the likelihood is the probability to die between last step unit time and current - step unit time, which is also the differences between probability to die before dh - and probability to die before dh-stepm . + /* i.e. if s2 is a death state and if the date of death is known + then the contribution to the likelihood is the probability to + die between last step unit time and current step unit time, + which is also equal to probability to die before dh + minus probability to die before dh-stepm . In version up to 0.92 likelihood was computed as if date of death was unknown. Death was treated as any other health state: the date of the interview describes the actual state @@ -1267,7 +1342,28 @@ double func( double *x) lower mortality. */ lli=log(out[s1][s2] - savm[s1][s2]); - }else{ + + + } else if (s2==-2) { + for (j=1,survp=0. ; j<=nlstate; j++) + survp += out[s1][j]; + lli= survp; + } + + else if (s2==-4) { + for (j=3,survp=0. ; j<=nlstate; j++) + survp += out[s1][j]; + lli= survp; + } + + else if (s2==-5) { + for (j=1,survp=0. ; j<=2; j++) + survp += out[s1][j]; + lli= survp; + } + + + else{ lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ /* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ } @@ -1831,7 +1927,7 @@ void lubksb(double **a, int n, int *indx } /************ Frequencies ********************/ -void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint) +void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[]) { /* Some frequencies */ int i, m, jk, k1,i1, j1, bool, z1,z2,j; @@ -1851,7 +1947,7 @@ void freqsummary(char fileres[], int ia fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); exit(0); } - freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3); + freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3); j1=0; j=cptcoveff; @@ -1864,8 +1960,8 @@ void freqsummary(char fileres[], int ia j1++; /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); scanf("%d", i);*/ - for (i=-1; i<=nlstate+ndeath; i++) - for (jk=-1; jk<=nlstate+ndeath; jk++) + for (i=-5; i<=nlstate+ndeath; i++) + for (jk=-5; jk<=nlstate+ndeath; jk++) for(m=iagemin; m <= iagemax+3; m++) freq[i][jk][m]=0; @@ -1904,7 +2000,7 @@ void freqsummary(char fileres[], int ia } /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ - +fprintf(ficresp, "#Local time at start: %s", strstart); if (cptcovn>0) { fprintf(ficresp, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); @@ -1990,7 +2086,7 @@ void freqsummary(char fileres[], int ia dateintmean=dateintsum/k2cpt; fclose(ficresp); - free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3); + free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3); free_vector(pp,1,nlstate); free_matrix(prop,1,nlstate,iagemin, iagemax+3); /* End of Freq */ @@ -2099,7 +2195,7 @@ void concatwav(int wav[], int **dh, int mi=0; m=firstpass; while(s[m][i] <= nlstate){ - if(s[m][i]>=1) + if(s[m][i]>=1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5) mw[++mi][i]=m; if(m >=lastpass) break; @@ -2117,11 +2213,11 @@ void concatwav(int wav[], int **dh, int if(mi==0){ nbwarn++; if(first==0){ - printf("Warning! None valid information for:%ld line=%d (skipped) and may be others, see log file\n",num[i],i); + printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i); first=1; } if(first==1){ - fprintf(ficlog,"Warning! None valid information for:%ld line=%d (skipped)\n",num[i],i); + fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i); } } /* end mi==0 */ } /* End individuals */ @@ -2139,13 +2235,19 @@ void concatwav(int wav[], int **dh, int nberr++; printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); j=1; /* Temporary Dangerous patch */ - printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n You MUST fix the contradiction between dates.\n",stepm); + printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); - fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n You MUST fix the contradiction between dates.\n",stepm); + fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); } k=k+1; - if (j >= jmax) jmax=j; - if (j <= jmin) jmin=j; + if (j >= jmax){ + jmax=j; + ijmax=i; + } + if (j <= jmin){ + jmin=j; + ijmin=i; + } sum=sum+j; /*if (j<0) printf("j=%d num=%d \n",j,i);*/ /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ @@ -2153,10 +2255,17 @@ void concatwav(int wav[], int **dh, int } else{ j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); - /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ +/* if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */ + k=k+1; - if (j >= jmax) jmax=j; - else if (j <= jmin)jmin=j; + if (j >= jmax) { + jmax=j; + ijmax=i; + } + else if (j <= jmin){ + jmin=j; + ijmin=i; + } /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/ if(j<0){ @@ -2199,8 +2308,8 @@ void concatwav(int wav[], int **dh, int } /* end wave */ } jmean=sum/k; - printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean); - fprintf(ficlog,"Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean); + printf("Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, num[ijmin], jmax, num[ijmax], jmean); + fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean); } /*********** Tricode ****************************/ @@ -2246,7 +2355,7 @@ void tricode(int *Tvar, int **nbcode, in for (k=0; k< maxncov; k++) Ndum[k]=0; for (i=1; i<=ncovmodel-2; i++) { - /* Listing of all covariables in staement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ + /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ ij=Tvar[i]; Ndum[ij]++; } @@ -2264,7 +2373,7 @@ void tricode(int *Tvar, int **nbcode, in /*********** Health Expectancies ****************/ -void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm,double delti[],double **matcov ) +void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm,double delti[],double **matcov,char strstart[] ) { /* Health expectancies */ @@ -2282,6 +2391,7 @@ void evsij(char fileres[], double ***eij dnewm=matrix(1,nlstate*nlstate,1,npar); doldm=matrix(1,nlstate*nlstate,1,nlstate*nlstate); + fprintf(ficreseij,"# Local time at start: %s", strstart); fprintf(ficreseij,"# Health expectancies\n"); fprintf(ficreseij,"# Age"); for(i=1; i<=nlstate;i++) @@ -2436,7 +2546,7 @@ void evsij(char fileres[], double ***eij } /************ Variance ******************/ -void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav) +void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[]) { /* Variance of health expectancies */ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ @@ -2487,7 +2597,9 @@ void varevsij(char optionfilefiname[], d fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev); } printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); + fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); + 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,"# Age cov=%-d",ij); for(j=nlstate+1; j<=(nlstate+ndeath);j++){ @@ -2497,11 +2609,12 @@ void varevsij(char optionfilefiname[], d } fprintf(ficresprobmorprev,"\n"); fprintf(ficgp,"\n# Routine varevsij"); + /* fprintf(fichtm, "#Local time at start: %s", strstart);*/ fprintf(fichtm,"\n

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

    • \n"); fprintf(fichtm,"\n
    %s
    \n",digitp); /* } */ varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); - + fprintf(ficresvij, "#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,"# Age"); for(i=1; i<=nlstate;i++) @@ -2734,7 +2847,7 @@ void varevsij(char optionfilefiname[], d } /* end varevsij */ /************ Variance of prevlim ******************/ -void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij) +void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, char strstart[]) { /* Variance of prevalence limit */ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ @@ -2747,7 +2860,7 @@ void varprevlim(char fileres[], double * double **gradg, **trgradg; double age,agelim; int theta; - + fprintf(ficresvpl, "#Local time at start: %s", strstart); fprintf(ficresvpl,"# Standard deviation of stable prevalences \n"); fprintf(ficresvpl,"# Age"); for(i=1; i<=nlstate;i++) @@ -2817,7 +2930,7 @@ void varprevlim(char fileres[], double * } /************ Variance of one-step probabilities ******************/ -void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax) +void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax, char strstart[]) { int i, j=0, i1, k1, l1, t, tj; int k2, l2, j1, z1; @@ -2862,11 +2975,13 @@ void varprob(char optionfilefiname[], do fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); - + fprintf(ficresprob, "#Local time at start: %s", strstart); fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n"); fprintf(ficresprob,"# Age"); + fprintf(ficresprobcov, "#Local time at start: %s", strstart); fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); fprintf(ficresprobcov,"# Age"); + fprintf(ficresprobcor, "#Local time at start: %s", strstart); fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); fprintf(ficresprobcov,"# Age"); @@ -3137,7 +3252,10 @@ void printinghtml(char fileres[], char t double jprev2, double mprev2,double anprev2){ int jj1, k1, i1, cpt; - fprintf(fichtm,"

    • Result files (first order: no variance)

      \n \ + fprintf(fichtm,""); + fprintf(fichtm,"

      • Result files (first order: no variance)

        \n \ - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): %s
        \n ", jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p")); fprintf(fichtm,"\ @@ -3188,7 +3306,7 @@ fprintf(fichtm," \n
        • Graphs fprintf(fichtm,"\ -\n

        • Result files (second order: variances)

          \n\ +\n

        • Result files (second order: variances)

          \n\ - Parameter file with estimated parameters and covariance matrix: %s
          \n", rfileres,rfileres); fprintf(fichtm," - Variance of one-step probabilities: %s
          \n", @@ -3902,28 +4020,28 @@ 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; + /* sump=sump+1;*/ num=num+1; } - /* for (i=1; i<=imx; i++) + /* for (i=0; i<=imx; i++) if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/ - for (i=0;i<=imx-1 ; i++) + for (i=1;i<=imx ; i++) { - if (cens[i]==1 & wav[i]>1) - A=-x[1]/(x[2])* - (exp(x[2]/YEARM*(agecens[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12))); + if (cens[i] == 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]/YEARM*(agedc[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12))) - +log(x[1]/YEARM)+x[2]/YEARM*(agedc[i]*12-agegomp*12)+log(YEARM); + if (cens[i] == 0 && wav[i]>1) + A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))) + +log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM); - if (wav[i]>1 & agecens[i]>15) { + /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */ + if (wav[i] > 1 ) { /* ??? */ L=L+A*weight[i]; /* printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/ } @@ -3937,8 +4055,8 @@ double gompertz(double x[]) /******************* Printing html file ***********/ void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \ int lastpass, int stepm, int weightopt, char model[],\ - int imx, double p[],double **matcov){ - int i; + int imx, double p[],double **matcov,double agemortsup){ + int i,k; fprintf(fichtm,"

          • Result files

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

            ",p[1],p[2],agegomp); @@ -3946,6 +4064,15 @@ void printinghtmlmort(char fileres[], ch fprintf(fichtm," p[%d] = %lf [%f ; %f]
            \n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); fprintf(fichtm,"

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

          • Life table

            \n
            "); + + fprintf(fichtm,"\nAge lx qx d(x,x+1) Lx Tx e
            "); + + for (k=agegomp;k<(agemortsup-2);k++) + fprintf(fichtm,"%d %.0lf %lf %.0lf %.0lf %.0lf %lf
            \n",k,lsurv[k],p[1]*exp(p[2]*(k-agegomp)),(p[1]*exp(p[2]*(k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]); + + fflush(fichtm); } @@ -3983,14 +4110,17 @@ int main(int argc, char *argv[]) { int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav); int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod; + int linei, month, year,iout; int jj, ll, li, lj, lk, imk; int numlinepar=0; /* Current linenumber of parameter file */ int itimes; int NDIM=2; char ca[32], cb[32], cc[32]; + char dummy[]=" "; /* FILE *fichtm; *//* Html File */ /* FILE *ficgp;*/ /*Gnuplot File */ + struct stat info; double agedeb, agefin,hf; double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; @@ -4003,7 +4133,7 @@ int main(int argc, char *argv[]) int *indx; char line[MAXLINE], linepar[MAXLINE]; char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; - char pathr[MAXLINE]; + char pathr[MAXLINE], pathimach[MAXLINE]; int firstobs=1, lastobs=10; int sdeb, sfin; /* Status at beginning and end */ int c, h , cpt,l; @@ -4013,6 +4143,8 @@ int main(int argc, char *argv[]) int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ int mobilav=0,popforecast=0; int hstepm, nhstepm; + int agemortsup; + float sumlpop=0.; double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000; double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000; @@ -4089,8 +4221,14 @@ int main(int argc, char *argv[]) printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/ /* cutv(path,optionfile,pathtot,'\\');*/ + /* Split argv[0], imach program to get pathimach */ + printf("\nargv[0]=%s argv[1]=%s, \n",argv[0],argv[1]); + split(argv[0],pathimach,optionfile,optionfilext,optionfilefiname); + printf("\nargv[0]=%s pathimach=%s, \noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",argv[0],pathimach,optionfile,optionfilext,optionfilefiname); + /* strcpy(pathimach,argv[0]); */ + /* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */ split(pathtot,path,optionfile,optionfilext,optionfilefiname); - printf("pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); + printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); chdir(path); strcpy(command,"mkdir "); strcat(command,optionfilefiname); @@ -4115,12 +4253,12 @@ int main(int argc, char *argv[]) } fprintf(ficlog,"Log filename:%s\n",filelog); fprintf(ficlog,"\n%s\n%s",version,fullversion); - fprintf(ficlog,"\nEnter the parameter file name: "); - fprintf(ficlog,"pathtot=%s\n\ + fprintf(ficlog,"\nEnter the parameter file name: \n"); + fprintf(ficlog,"pathimach=%s\npathtot=%s\n\ path=%s \n\ optionfile=%s\n\ optionfilext=%s\n\ - optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); + optionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname); printf("Local time (at start):%s",strstart); fprintf(ficlog,"Local time (at start): %s",strstart); @@ -4257,7 +4395,6 @@ int main(int argc, char *argv[]) } fflush(ficlog); - p=param[1][1]; /* Reads comments: lines beginning with '#' */ @@ -4381,40 +4518,116 @@ int main(int argc, char *argv[]) ncodemax=ivector(1,8); i=1; - while (fgets(line, MAXLINE, fic) != NULL) { - if ((i >= firstobs) && (i <=lastobs)) { - - for (j=maxwav;j>=1;j--){ - cutv(stra, strb,line,' '); s[j][i]=atoi(strb); - strcpy(line,stra); - cutv(stra, strb,line,'/'); anint[j][i]=(double)(atoi(strb)); strcpy(line,stra); - cutv(stra, strb,line,' '); mint[j][i]=(double)(atoi(strb)); strcpy(line,stra); + linei=0; + while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) { + linei=linei+1; + for(j=strlen(line); j>=0;j--){ /* Untabifies line */ + if(line[j] == '\t') + line[j] = ' '; + } + for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10)||(line[j]==13);j--){ + ; + }; + line[j+1]=0; /* Trims blanks at end of line */ + if(line[0]=='#'){ + fprintf(ficlog,"Comment line\n%s\n",line); + printf("Comment line\n%s\n",line); + continue; + } + + for (j=maxwav;j>=1;j--){ + cutv(stra, strb,line,' '); + errno=0; + lval=strtol(strb,&endptr,10); + /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/ + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%d' at line number %d %s for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav); + exit(1); } - - cutv(stra, strb,line,'/'); andc[i]=(double)(atoi(strb)); strcpy(line,stra); - cutv(stra, strb,line,' '); moisdc[i]=(double)(atoi(strb)); strcpy(line,stra); - - cutv(stra, strb,line,'/'); annais[i]=(double)(atoi(strb)); strcpy(line,stra); - cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra); - - cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra); - for (j=ncovcol;j>=1;j--){ - cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra); - } - lstra=strlen(stra); - if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */ - stratrunc = &(stra[lstra-9]); - num[i]=atol(stratrunc); + s[j][i]=lval; + + strcpy(line,stra); + cutv(stra, strb,line,' '); + if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ + } + else if(iout=sscanf(strb,"%s.") != 0){ + month=99; + year=9999; + }else{ + printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j); + exit(1); + } + anint[j][i]= (double) year; + mint[j][i]= (double)month; + strcpy(line,stra); + } /* ENd Waves */ + + cutv(stra, strb,line,' '); + if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ + } + else if(iout=sscanf(strb,"%s.",dummy) != 0){ + month=99; + year=9999; + }else{ + printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line); + exit(1); + } + andc[i]=(double) year; + moisdc[i]=(double) month; + strcpy(line,stra); + + cutv(stra, strb,line,' '); + if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ + } + else if(iout=sscanf(strb,"%s.") != 0){ + month=99; + year=9999; + }else{ + printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j); + exit(1); + } + annais[i]=(double)(year); + moisnais[i]=(double)(month); + strcpy(line,stra); + + cutv(stra, strb,line,' '); + errno=0; + lval=strtol(strb,&endptr,10); + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a weight. Exiting.\n",lval, i,line,linei); + exit(1); + } + weight[i]=(double)(lval); + strcpy(line,stra); + + for (j=ncovcol;j>=1;j--){ + cutv(stra, strb,line,' '); + errno=0; + lval=strtol(strb,&endptr,10); + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a covar (meaning 0 for the reference or 1). Exiting.\n",lval, linei,i, line); + exit(1); + } + if(lval <-1 || lval >1){ + printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a value of the %d covar (meaning 0 for the reference or 1. IMaCh does not build design variables, do it your self). Exiting.\n",lval,linei, i,line,j); + exit(1); } - else - num[i]=atol(stra); - - /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){ - printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/ - - i=i+1; - } - } + covar[j][i]=(double)(lval); + strcpy(line,stra); + } + lstra=strlen(stra); + + if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */ + stratrunc = &(stra[lstra-9]); + num[i]=atol(stratrunc); + } + else + num[i]=atol(stra); + /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){ + printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/ + + i=i+1; + } /* End loop reading data */ /* printf("ii=%d", ij); scanf("%d",i);*/ imx=i-1; /* Number of individuals */ @@ -4428,12 +4641,12 @@ int main(int argc, char *argv[]) if (s[4][i]==9) s[4][i]=-1; printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i]));}*/ - for (i=1; i<=imx; i++) + /* for (i=1; i<=imx; i++) */ /*if ((s[3][i]==3) || (s[4][i]==3)) weight[i]=0.08; else weight[i]=1;*/ - /* Calculation of the number of parameter from char model*/ + /* Calculation of the number of parameters from char model */ Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */ Tprod=ivector(1,15); Tvaraff=ivector(1,15); @@ -4546,7 +4759,7 @@ int main(int argc, char *argv[]) for (i=1; i<=imx; i++) { agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); for(m=firstpass; (m<= lastpass); m++){ - if(s[m][i] >0){ + if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){ if (s[m][i] >= nlstate+1) { if(agedc[i]>0) if((int)moisdc[i]!=99 && (int)andc[i]!=9999) @@ -4562,8 +4775,7 @@ int main(int argc, char *argv[]) } } else if(s[m][i] !=9){ /* Standard case, age in fractional - years but with the precision of a - month */ + years but with the precision of a month */ agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]); if((int)mint[m][i]==99 || (int)anint[m][i]==9999) agev[m][i]=1; @@ -4720,7 +4932,7 @@ Title=%s
            Datafile=%s Firstpass=%d La /* Calculates basic frequencies. Computes observed prevalence at single age and prints on file fileres'p'. */ - freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart); fprintf(fichtm,"\n"); fprintf(fichtm,"
            Total number of observations=%d
            \n\ @@ -4739,6 +4951,7 @@ Interval (in months) between two waves: p=param[1][1]; /* *(*(*(param +1)+1)+0) */ globpr=0; /* To get the number ipmx of contributions and the sum of weights*/ + if (mle==-3){ ximort=matrix(1,NDIM,1,NDIM); cens=ivector(1,n); @@ -4748,9 +4961,9 @@ Interval (in months) between two waves: for (i=1; i<=imx; i++){ dcwave[i]=-1; - for (j=1; j<=lastpass; j++) - if (s[j][i]>nlstate) { - dcwave[i]=j; + for (m=firstpass; m<=lastpass; m++) + if (s[m][i]>nlstate) { + dcwave[i]=m; /* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/ break; } @@ -4759,12 +4972,16 @@ Interval (in months) between two waves: 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 ; + j=wav[i]; + agecens[i]=1.; + + if (ageexmed[i]> 1 && wav[i] > 0){ + agecens[i]=agev[mw[j][i]][i]; + cens[i]= 1; + }else if (ageexmed[i]< 1) + cens[i]= -1; + if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass) + cens[i]=0 ; } else cens[i]=-1; } @@ -4773,29 +4990,29 @@ Interval (in months) between two waves: for (j=1;j<=NDIM;j++) ximort[i][j]=(i == j ? 1.0 : 0.0); } - - p[1]=0.1; p[2]=0.1; + + p[1]=0.0268; p[NDIM]=0.083; /*printf("%lf %lf", p[1], p[2]);*/ - printf("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"); - + printf("Powell\n"); fprintf(ficlog,"Powell\n"); + strcpy(filerespow,"pow-mort"); + strcat(filerespow,fileres); + if((ficrespow=fopen(filerespow,"w"))==NULL) { + printf("Problem with resultfile: %s\n", filerespow); + fprintf(ficlog,"Problem with resultfile: %s\n", filerespow); + } + fprintf(ficrespow,"# Powell\n# iter -2*LL"); + /* for (i=1;i<=nlstate;i++) + for(j=1;j<=nlstate+ndeath;j++) + if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); + */ + fprintf(ficrespow,"\n"); + powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); fclose(ficrespow); - hesscov(matcov, p, NDIM,delti, 1e-4, gompertz); + hesscov(matcov, p, NDIM, delti, 1e-4, gompertz); for(i=1; i <=NDIM; i++) for(j=i+1;j<=NDIM;j++) @@ -4812,14 +5029,49 @@ Interval (in months) between two waves: printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp); for (i=1;i<=NDIM;i++) printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); + + lsurv=vector(1,AGESUP); + lpop=vector(1,AGESUP); + tpop=vector(1,AGESUP); + lsurv[agegomp]=100000; + + for (k=agegomp;k<=AGESUP;k++) { + agemortsup=k; + if (p[1]*exp(p[2]*(k-agegomp))>1) break; + } + + for (k=agegomp;k=1 */ likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ @@ -5077,6 +5329,7 @@ Interval (in months) between two waves: } printf("Computing stable prevalence: result on file '%s' \n", filerespl); fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl); + fprintf(ficrespl, "#Local time at start: %s", strstart); fprintf(ficrespl,"#Stable prevalence \n"); fprintf(ficrespl,"#Age "); for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); @@ -5137,7 +5390,7 @@ Interval (in months) between two waves: hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ /* hstepm=1; aff par mois*/ - + fprintf(ficrespij, "#Local time at start: %s", strstart); fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x "); for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ @@ -5174,7 +5427,7 @@ Interval (in months) between two waves: } } - varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax); + varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart); fclose(ficrespij); @@ -5263,16 +5516,16 @@ Interval (in months) between two waves: 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); + evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart); vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); oldm=oldms;savm=savms; - varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav); + varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav, strstart); if(popbased==1){ - varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav); + varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav, strstart); } - + 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"); @@ -5345,7 +5598,7 @@ Interval (in months) between two waves: varpl=matrix(1,nlstate,(int) bage, (int) fage); oldm=oldms;savm=savms; - varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k); + 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); } } @@ -5394,7 +5647,7 @@ Interval (in months) between two waves: tm = *localtime(&end_time.tv_sec); tmg = *gmtime(&end_time.tv_sec); strcpy(strtend,asctime(&tm)); - printf("Local time at start %s\nLocaltime at end %s",strstart, strtend); + printf("Local time at start %s\nLocal time at end %s",strstart, strtend); fprintf(ficlog,"Local time at start %s\nLocal time at end %s\n",strstart, strtend); printf("Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout)); @@ -5411,12 +5664,32 @@ Interval (in months) between two waves: /*------ End -----------*/ chdir(path); - strcpy(plotcmd,GNUPLOTPROGRAM); - strcat(plotcmd," "); - strcat(plotcmd,optionfilegnuplot); - printf("Starting graphs with: %s",plotcmd);fflush(stdout); + /*strcat(plotcmd,CHARSEPARATOR);*/ + sprintf(plotcmd,"gnuplot"); +#ifndef UNIX + sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach); +#endif + if(!stat(plotcmd,&info)){ + printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout); + if(!stat(getenv("GNUPLOTBIN"),&info)){ + printf("Error gnuplot program not found: %s Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout); + }else + strcpy(pplotcmd,plotcmd); +#ifdef UNIX + strcpy(plotcmd,GNUPLOTPROGRAM); + if(!stat(plotcmd,&info)){ + printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout); + }else + strcpy(pplotcmd,plotcmd); +#endif + }else + strcpy(pplotcmd,plotcmd); + + sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot); + printf("Starting graphs with: %s\n",plotcmd);fflush(stdout); + if((outcmd=system(plotcmd)) != 0){ - printf(" Problem with gnuplot\n"); + printf("\n Problem with gnuplot\n"); } printf(" Wait..."); while (z[0] != 'q') { @@ -5424,7 +5697,10 @@ Interval (in months) between two waves: printf("\nType e to edit output files, g to graph again and q for exiting: "); scanf("%s",z); /* if (z[0] == 'c') system("./imach"); */ - if (z[0] == 'e') system(optionfilehtm); + if (z[0] == 'e') { + printf("Starting browser with: %s",optionfilehtm);fflush(stdout); + system(optionfilehtm); + } else if (z[0] == 'g') system(plotcmd); else if (z[0] == 'q') exit(0); }