--- imach/src/imach.c 2006/01/19 18:05:42 1.108 +++ imach/src/imach.c 2006/03/20 09:45:41 1.122 @@ -1,6 +1,74 @@ -/* $Id: imach.c,v 1.108 2006/01/19 18:05:42 lievre Exp $ +/* $Id: imach.c,v 1.122 2006/03/20 09:45:41 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.122 2006/03/20 09:45:41 brouard + (Module): Weights can have a decimal point as for + English (a comma might work with a correct LC_NUMERIC environment, + otherwise the weight is truncated). + Modification of warning when the covariates values are not 0 or + 1. + Version 0.98g + + Revision 1.121 2006/03/16 17:45:01 lievre + * imach.c (Module): Comments concerning covariates added + + * imach.c (Module): refinements in the computation of lli if + status=-2 in order to have more reliable computation if stepm is + not 1 month. Version 0.98f + + Revision 1.120 2006/03/16 15:10:38 lievre + (Module): refinements in the computation of lli if + status=-2 in order to have more reliable computation if stepm is + not 1 month. Version 0.98f + + Revision 1.119 2006/03/15 17:42:26 brouard + (Module): Bug if status = -2, the loglikelihood was + computed as likelihood omitting the logarithm. Version O.98e + + Revision 1.118 2006/03/14 18:20:07 brouard + (Module): varevsij Comments added explaining the second + table of variances if popbased=1 . + (Module): Covariances of eij, ekl added, graphs fixed, new html link. + (Module): Function pstamp added + (Module): Version 0.98d + + Revision 1.117 2006/03/14 17:16:22 brouard + (Module): varevsij Comments added explaining the second + table of variances if popbased=1 . + (Module): Covariances of eij, ekl added, graphs fixed, new html link. + (Module): Function pstamp added + (Module): Version 0.98d + + Revision 1.116 2006/03/06 10:29:27 brouard + (Module): Variance-covariance wrong links and + varian-covariance of ej. is needed (Saito). + + Revision 1.115 2006/02/27 12:17:45 brouard + (Module): One freematrix added in mlikeli! 0.98c + + Revision 1.114 2006/02/26 12:57:58 brouard + (Module): Some improvements in processing parameter + filename with strsep. + + Revision 1.113 2006/02/24 14:20:24 brouard + (Module): Memory leaks checks with valgrind and: + datafile was not closed, some imatrix were not freed and on matrix + allocation too. + + Revision 1.112 2006/01/30 09:55:26 brouard + (Module): Back to gnuplot.exe instead of wgnuplot.exe + + Revision 1.111 2006/01/25 20:38:18 brouard + (Module): Lots of cleaning and bugs added (Gompertz) + (Module): Comments can be added in data file. Missing date values + can be a simple dot '.'. + + Revision 1.110 2006/01/25 00:51:50 brouard + (Module): Lots of cleaning and bugs added (Gompertz) + + Revision 1.109 2006/01/24 19:37:15 brouard + (Module): Comments (lines starting with a #) are allowed in data. + Revision 1.108 2006/01/19 18:05:42 lievre Gnuplot problem appeared... To be fixed @@ -180,7 +248,7 @@ hPijx. Also this programme outputs the covariance matrix of the parameters but also - of the life expectancies. It also computes the stable prevalence. + of the life expectancies. It also computes the period (stable) prevalence. Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). Institut national d'études démographiques, Paris. @@ -209,7 +277,7 @@ begin-prev-date,... open gnuplot file open html file - stable prevalence + period (stable) prevalence for age prevalim() h Pij x variance of p varprob @@ -221,7 +289,7 @@ varevsij() if popbased==1 varevsij(,popbased) total life expectancies - Variance of stable prevalence + Variance of period (stable) prevalence end */ @@ -234,6 +302,7 @@ #include #include +#include #include #include #include @@ -277,11 +346,13 @@ extern int errno; #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.108 2006/01/19 18:05:42 lievre Exp $ */ +/* $Id: imach.c,v 1.122 2006/03/20 09:45:41 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.98a, January 2006, INED-EUROREVES "; -char fullversion[]="$Revision: 1.108 $ $Date: 2006/01/19 18:05:42 $"; +char version[]="Imach version 0.98g, March 2006, INED-EUROREVES-Institut de longevite "; +char fullversion[]="$Revision: 1.122 $ $Date: 2006/03/20 09:45:41 $"; +char strstart[80]; +char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ int nvar; int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; @@ -294,6 +365,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; @@ -304,7 +376,7 @@ int **bh; /* bh[mi][i] is the bias (+ or double jmean; /* Mean space between 2 waves */ double **oldm, **newm, **savm; /* Working pointers to matrices */ double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ -FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; +FILE *fic,*ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; FILE *ficlog, *ficrespow; int globpr; /* Global variable for printing or not */ double fretone; /* Only one call to likelihood */ @@ -318,13 +390,17 @@ FILE *ficresprobmorprev; FILE *fichtm, *fichtmcov; /* Html File */ FILE *ficreseij; char filerese[FILENAMELENGTH]; +FILE *ficresstdeij; +char fileresstde[FILENAMELENGTH]; +FILE *ficrescveij; +char filerescve[FILENAMELENGTH]; FILE *ficresvij; char fileresv[FILENAMELENGTH]; FILE *ficresvpl; char fileresvpl[FILENAMELENGTH]; char title[MAXLINE]; char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH]; -char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH]; +char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH]; char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH]; char command[FILENAMELENGTH]; int outcmd=0; @@ -346,6 +422,10 @@ long time_value; extern long time(); char strcurr[80], strfor[80]; +char *endptr; +long lval; +double dval; + #define NR_END 1 #define FREE_ARG char* #define FTOL 1.0e-10 @@ -1034,7 +1114,7 @@ void powell(double p[], double **xi, int } } -/**** Prevalence limit (stable prevalence) ****************/ +/**** Prevalence limit (stable or period prevalence) ****************/ double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij) { @@ -1327,23 +1407,23 @@ double func( double *x) } else if (s2==-2) { for (j=1,survp=0. ; j<=nlstate; j++) - survp += out[s1][j]; - lli= survp; + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + /*survp += out[s1][j]; */ + lli= log(survp); } - else if (s2==-4) { - for (j=3,survp=0. ; j<=nlstate; j++) - survp += out[s1][j]; - lli= survp; - } - - else if (s2==-5) { - for (j=1,survp=0. ; j<=2; j++) - survp += out[s1][j]; - lli= survp; - } - + else if (s2==-4) { + for (j=3,survp=0. ; j<=nlstate; j++) + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + lli= log(survp); + } + else if (s2==-5) { + for (j=1,survp=0. ; j<=2; j++) + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + lli= log(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 */ @@ -1538,7 +1618,11 @@ double funcone( double *x) */ if( s2 > nlstate && (mle <5) ){ /* Jackson */ lli=log(out[s1][s2] - savm[s1][s2]); - } else if (mle==1){ + } else if (s2==-2) { + for (j=1,survp=0. ; j<=nlstate; j++) + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + lli= log(survp); + }else if (mle==1){ lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ } else if(mle==2){ lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ @@ -1554,8 +1638,8 @@ double funcone( double *x) ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; /* printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ if(globpr){ - fprintf(ficresilk,"%9d %6d %1d %1d %1d %1d %3d %10.6f %6.4f\ - %10.6f %10.6f %10.6f ", \ + fprintf(ficresilk,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\ + %11.6f %11.6f %11.6f ", \ num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i], 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); for(k=1,llt=0.,l=0.; k<=nlstate; k++){ @@ -1640,6 +1724,7 @@ void mlikeli(FILE *ficres,double p[], in powell(p,xi,npar,ftol,&iter,&fret,func); + free_matrix(xi,1,npar,1,npar); fclose(ficrespow); printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p)); fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); @@ -1907,6 +1992,11 @@ void lubksb(double **a, int n, int *indx } } +void pstamp(FILE *fichier) +{ + fprintf(fichier,"# %s.%s\n#%s\n#%s\n# %s", optionfilefiname,optionfilext,version,fullversion,strstart); +} + /************ Frequencies ********************/ void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[]) { /* Some frequencies */ @@ -1916,7 +2006,6 @@ void freqsummary(char fileres[], int ia double ***freq; /* Frequencies */ double *pp, **prop; double pos,posprop, k2, dateintsum=0,k2cpt=0; - FILE *ficresp; char fileresp[FILENAMELENGTH]; pp=vector(1,nlstate); @@ -1981,7 +2070,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); + pstamp(ficresp); if (cptcovn>0) { fprintf(ficresp, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); @@ -2194,11 +2283,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 */ @@ -2221,8 +2310,14 @@ void concatwav(int wav[], int **dh, int 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);*/ @@ -2233,8 +2328,14 @@ void concatwav(int wav[], int **dh, int /* 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){ @@ -2277,8 +2378,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 ****************************/ @@ -2342,32 +2443,151 @@ void tricode(int *Tvar, int **nbcode, in /*********** Health Expectancies ****************/ -void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm,double delti[],double **matcov,char strstart[] ) +void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] ) + +{ + /* Health expectancies, no variances */ + int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2; + double age, agelim, hf; + double ***p3mat; + double eip; + + pstamp(ficreseij); + fprintf(ficreseij,"# (a) Life expectancies by health status at initial age and (b) health expectancies by health status at initial age\n"); + fprintf(ficreseij,"# Age"); + for(i=1; i<=nlstate;i++){ + for(j=1; j<=nlstate;j++){ + fprintf(ficreseij," e%1d%1d ",i,j); + } + fprintf(ficreseij," e%1d. ",i); + } + fprintf(ficreseij,"\n"); + + + if(estepm < stepm){ + printf ("Problem %d lower than %d\n",estepm, stepm); + } + else hstepm=estepm; + /* We compute the life expectancy from trapezoids spaced every estepm months + * This is mainly to measure the difference between two models: for example + * if stepm=24 months pijx are given only every 2 years and by summing them + * we are calculating an estimate of the Life Expectancy assuming a linear + * progression in between and thus overestimating or underestimating according + * to the curvature of the survival function. If, for the same date, we + * estimate the model with stepm=1 month, we can keep estepm to 24 months + * to compare the new estimate of Life expectancy with the same linear + * hypothesis. A more precise result, taking into account a more precise + * curvature will be obtained if estepm is as small as stepm. */ + + /* For example we decided to compute the life expectancy with the smallest unit */ + /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. + nhstepm is the number of hstepm from age to agelim + nstepm is the number of stepm from age to agelin. + Look at hpijx to understand the reason of that which relies in memory size + and note for a fixed period like estepm months */ + /* We decided (b) to get a life expectancy respecting the most precise curvature of the + survival function given by stepm (the optimization length). Unfortunately it + means that if the survival funtion is printed only each two years of age and if + you sum them up and add 1 year (area under the trapezoids) you won't get the same + results. So we changed our mind and took the option of the best precision. + */ + hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ + + agelim=AGESUP; + /* nhstepm age range expressed in number of stepm */ + nstepm=(int) rint((agelim-age)*YEARM/stepm); + /* Typically if 20 years nstepm = 20*12/6=40 stepm */ + /* if (stepm >= YEARM) hstepm=1;*/ + nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ + p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + + for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ + /* Computed by stepm unit matrices, product of hstepm matrices, stored + in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ + + hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); + + hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ + + printf("%d|",(int)age);fflush(stdout); + fprintf(ficlog,"%d|",(int)age);fflush(ficlog); + + /* Computing expectancies */ + for(i=1; i<=nlstate;i++) + for(j=1; j<=nlstate;j++) + for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ + eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf; + + /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ + + } + + fprintf(ficreseij,"%3.0f",age ); + for(i=1; i<=nlstate;i++){ + eip=0; + for(j=1; j<=nlstate;j++){ + eip +=eij[i][j][(int)age]; + fprintf(ficreseij,"%9.4f", eij[i][j][(int)age] ); + } + fprintf(ficreseij,"%9.4f", eip ); + } + fprintf(ficreseij,"\n"); + + } + free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + printf("\n"); + fprintf(ficlog,"\n"); + +} + +void cvevsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] ) { - /* Health expectancies */ - int i, j, nhstepm, hstepm, h, nstepm, k, cptj; + /* Covariances of health expectancies eij and of total life expectancies according + to initial status i, ei. . + */ + int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji; double age, agelim, hf; - double ***p3mat,***varhe; + double ***p3matp, ***p3matm, ***varhe; double **dnewm,**doldm; - double *xp; + double *xp, *xm; double **gp, **gm; double ***gradg, ***trgradg; int theta; + double eip, vip; + varhe=ma3x(1,nlstate*nlstate,1,nlstate*nlstate,(int) bage, (int) fage); xp=vector(1,npar); + xm=vector(1,npar); dnewm=matrix(1,nlstate*nlstate,1,npar); doldm=matrix(1,nlstate*nlstate,1,nlstate*nlstate); - fprintf(ficreseij,"# Local time at start: %s", strstart); - fprintf(ficreseij,"# Health expectancies\n"); - fprintf(ficreseij,"# Age"); - for(i=1; i<=nlstate;i++) + pstamp(ficresstdeij); + fprintf(ficresstdeij,"# Health expectancies with standard errors\n"); + fprintf(ficresstdeij,"# Age"); + for(i=1; i<=nlstate;i++){ for(j=1; j<=nlstate;j++) - fprintf(ficreseij," %1d-%1d (SE)",i,j); - fprintf(ficreseij,"\n"); + fprintf(ficresstdeij," e%1d%1d (SE)",i,j); + fprintf(ficresstdeij," e%1d. ",i); + } + fprintf(ficresstdeij,"\n"); + pstamp(ficrescveij); + fprintf(ficrescveij,"# Subdiagonal matrix of covariances of health expectancies by age: cov(eij,ekl)\n"); + fprintf(ficrescveij,"# Age"); + for(i=1; i<=nlstate;i++) + for(j=1; j<=nlstate;j++){ + cptj= (j-1)*nlstate+i; + for(i2=1; i2<=nlstate;i2++) + for(j2=1; j2<=nlstate;j2++){ + cptj2= (j2-1)*nlstate+i2; + if(cptj2 <= cptj) + fprintf(ficrescveij," %1d%1d,%1d%1d",i,j,i2,j2); + } + } + fprintf(ficrescveij,"\n"); + if(estepm < stepm){ printf ("Problem %d lower than %d\n",estepm, stepm); } @@ -2397,77 +2617,64 @@ void evsij(char fileres[], double ***eij */ hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ + /* If stepm=6 months */ + /* nhstepm age range expressed in number of stepm */ agelim=AGESUP; - for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ - /* nhstepm age range expressed in number of stepm */ - nstepm=(int) rint((agelim-age)*YEARM/stepm); - /* Typically if 20 years nstepm = 20*12/6=40 stepm */ - /* if (stepm >= YEARM) hstepm=1;*/ - nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ - p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - gradg=ma3x(0,nhstepm,1,npar,1,nlstate*nlstate); - gp=matrix(0,nhstepm,1,nlstate*nlstate); - gm=matrix(0,nhstepm,1,nlstate*nlstate); + nstepm=(int) rint((agelim-age)*YEARM/stepm); + /* Typically if 20 years nstepm = 20*12/6=40 stepm */ + /* if (stepm >= YEARM) hstepm=1;*/ + nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ + + p3matp=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + p3matm=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + gradg=ma3x(0,nhstepm,1,npar,1,nlstate*nlstate); + trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar); + gp=matrix(0,nhstepm,1,nlstate*nlstate); + gm=matrix(0,nhstepm,1,nlstate*nlstate); + + for (age=bage; age<=fage; age ++){ /* Computed by stepm unit matrices, product of hstepm matrices, stored in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ - hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij); - hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ /* Computing Variances of health expectancies */ - - for(theta=1; theta <=npar; theta++){ + /* Gradient is computed with plus gp and minus gm. Code is duplicated in order to + decrease memory allocation */ + for(theta=1; theta <=npar; theta++){ for(i=1; i<=npar; i++){ xp[i] = x[i] + (i==theta ?delti[theta]:0); + xm[i] = x[i] - (i==theta ?delti[theta]:0); } - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); + hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij); + hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij); - cptj=0; for(j=1; j<= nlstate; j++){ for(i=1; i<=nlstate; i++){ - cptj=cptj+1; - for(h=0, gp[h][cptj]=0.; h<=nhstepm-1; h++){ - gp[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.; + for(h=0; h<=nhstepm-1; h++){ + gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.; + gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.; } } } - - for(i=1; i<=npar; i++) - xp[i] = x[i] - (i==theta ?delti[theta]:0); - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); - - cptj=0; - for(j=1; j<= nlstate; j++){ - for(i=1;i<=nlstate;i++){ - cptj=cptj+1; - for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){ - - gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.; - } - } - } - for(j=1; j<= nlstate*nlstate; j++) + for(ij=1; ij<= nlstate*nlstate; ij++) for(h=0; h<=nhstepm-1; h++){ - gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; + gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta]; } - } - -/* End theta */ - - trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar); - - for(h=0; h<=nhstepm-1; h++) + }/* End theta */ + + + for(h=0; h<=nhstepm-1; h++) for(j=1; j<=nlstate*nlstate;j++) for(theta=1; theta <=npar; theta++) trgradg[h][j][theta]=gradg[h][theta][j]; - + - for(i=1;i<=nlstate*nlstate;i++) - for(j=1;j<=nlstate*nlstate;j++) - varhe[i][j][(int)age] =0.; + for(ij=1;ij<=nlstate*nlstate;ij++) + for(ji=1;ji<=nlstate*nlstate;ji++) + varhe[ij][ji][(int)age] =0.; printf("%d|",(int)age);fflush(stdout); fprintf(ficlog,"%d|",(int)age);fflush(ficlog); @@ -2475,39 +2682,60 @@ void evsij(char fileres[], double ***eij for(k=0;k<=nhstepm-1;k++){ matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov); matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]); - for(i=1;i<=nlstate*nlstate;i++) - for(j=1;j<=nlstate*nlstate;j++) - varhe[i][j][(int)age] += doldm[i][j]*hf*hf; + for(ij=1;ij<=nlstate*nlstate;ij++) + for(ji=1;ji<=nlstate*nlstate;ji++) + varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf; } } /* Computing expectancies */ + hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ - eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf; + eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[i][j][h+1])/2.0*hf; -/* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ + /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ } - fprintf(ficreseij,"%3.0f",age ); - cptj=0; + fprintf(ficresstdeij,"%3.0f",age ); + for(i=1; i<=nlstate;i++){ + eip=0.; + vip=0.; + for(j=1; j<=nlstate;j++){ + eip += eij[i][j][(int)age]; + for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */ + vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age]; + fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) ); + } + fprintf(ficresstdeij," %9.4f (%.4f)", eip, sqrt(vip)); + } + fprintf(ficresstdeij,"\n"); + + fprintf(ficrescveij,"%3.0f",age ); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++){ - cptj++; - fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) ); + cptj= (j-1)*nlstate+i; + for(i2=1; i2<=nlstate;i2++) + for(j2=1; j2<=nlstate;j2++){ + cptj2= (j2-1)*nlstate+i2; + if(cptj2 <= cptj) + fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]); + } } - fprintf(ficreseij,"\n"); + fprintf(ficrescveij,"\n"); - free_matrix(gm,0,nhstepm,1,nlstate*nlstate); - free_matrix(gp,0,nhstepm,1,nlstate*nlstate); - free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*nlstate); - free_ma3x(trgradg,0,nhstepm,1,nlstate*nlstate,1,npar); - free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); } + free_matrix(gm,0,nhstepm,1,nlstate*nlstate); + free_matrix(gp,0,nhstepm,1,nlstate*nlstate); + free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*nlstate); + free_ma3x(trgradg,0,nhstepm,1,nlstate*nlstate,1,npar); + free_ma3x(p3matm,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + free_ma3x(p3matp,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); printf("\n"); fprintf(ficlog,"\n"); + free_vector(xm,1,npar); free_vector(xp,1,npar); free_matrix(dnewm,1,nlstate*nlstate,1,npar); free_matrix(doldm,1,nlstate*nlstate,1,nlstate*nlstate); @@ -2568,7 +2796,7 @@ void varevsij(char optionfilefiname[], d 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); + pstamp(ficresprobmorprev); fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm); fprintf(ficresprobmorprev,"# Age cov=%-d",ij); for(j=nlstate+1; j<=(nlstate+ndeath);j++){ @@ -2583,12 +2811,16 @@ void varevsij(char optionfilefiname[], d 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"); + pstamp(ficresvij); + fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are "); + if(popbased==1) + fprintf(ficresvij,"the age specific prevalence observed in the population i.e cross-sectionally\n in each health state (popbased=1)"); + else + fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n"); fprintf(ficresvij,"# Age"); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) - fprintf(ficresvij," Cov(e%1d, e%1d)",i,j); + fprintf(ficresvij," Cov(e.%1d, e.%1d)",i,j); fprintf(ficresvij,"\n"); xp=vector(1,npar); @@ -2829,8 +3061,9 @@ 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"); + + pstamp(ficresvpl); + fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n"); fprintf(ficresvpl,"# Age"); for(i=1; i<=nlstate;i++) fprintf(ficresvpl," %1d-%1d",i,i); @@ -2944,15 +3177,15 @@ void varprob(char optionfilefiname[], do fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); - fprintf(ficresprob, "#Local time at start: %s", strstart); + pstamp(ficresprob); fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n"); fprintf(ficresprob,"# Age"); - fprintf(ficresprobcov, "#Local time at start: %s", strstart); + pstamp(ficresprobcov); fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); fprintf(ficresprobcov,"# Age"); - fprintf(ficresprobcor, "#Local time at start: %s", strstart); + pstamp(ficresprobcor); fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); - fprintf(ficresprobcov,"# Age"); + fprintf(ficresprobcor,"# Age"); for(i=1; i<=nlstate;i++) @@ -3203,6 +3436,8 @@ To be simple, these graphs help to under } free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); + free_matrix(doldm,1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); + free_matrix(dnewm,1,(nlstate)*(nlstate+ndeath),1,npar); free_vector(xp,1,npar); fclose(ficresprob); fclose(ficresprobcov); @@ -3231,13 +3466,14 @@ void printinghtml(char fileres[], char t - Estimated transition probabilities over %d (stepm) months: %s
\n ", stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij")); fprintf(fichtm,"\ - - Stable prevalence in each health state: %s
\n", + - Period (stable) prevalence in each health state: %s
\n", subdirf2(fileres,"pl"),subdirf2(fileres,"pl")); fprintf(fichtm,"\ - - Life expectancies by age and initial health status (estepm=%2d months): \ + - (a) Life expectancies by health status at initial age, (b) health expectancies by health status at initial age: ei., eij . If one or more covariate are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ %s
\n", estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e")); + fprintf(fichtm," \n
  • Graphs

    • "); m=cptcoveff; @@ -3254,20 +3490,20 @@ fprintf(fichtm," \n

    • Graphs fprintf(fichtm," ************\n
      "); } /* Pij */ - fprintf(fichtm,"
      - Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s%d1.png
      \ -",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); + fprintf(fichtm,"
      - Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s%d1.png
      \ +",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); /* Quasi-incidences */ fprintf(fichtm,"
      - Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\ - before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: %s%d2.png
      \ -",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); - /* Stable prevalence in each health state */ + before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: %s%d2.png
      \ +",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); + /* Period (stable) prevalence in each health state */ for(cpt=1; cpt- Stable prevalence in each health state : p%s%d%d.png
      \ -",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1); + fprintf(fichtm,"
      - Period (stable) prevalence in each health state : %s%d%d.png
      \ +",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1); } for(cpt=1; cpt<=nlstate;cpt++) { - fprintf(fichtm,"\n
      - Health life expectancies by age and initial health state (%d): %s%d%d.png
      \ -",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1); + fprintf(fichtm,"\n
      - Life expectancy by health state (%d) at initial age and its decomposition into health expectancies : %s%d%d.png
      \ +",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1); } } /* end i1 */ }/* End k1 */ @@ -3288,13 +3524,21 @@ fprintf(fichtm," \n
      • Graphs - Correlation matrix of one-step probabilities: %s
        \n", subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor")); fprintf(fichtm,"\ - - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): %s
        \n", + - Variances and covariances of health expectancies by age and initial health status (cov(eij,ekl)(estepm=%2d months): \ + %s
        \n
        • ", + estepm,subdirf2(fileres,"cve"),subdirf2(fileres,"cve")); + fprintf(fichtm,"\ + - (a) Health expectancies by health status at initial age (eij) and standard errors (in parentheses) (b) life expectancies and standard errors (ei.=ei1+ei2+...)(estepm=%2d months): \ + %s
        \n", + estepm,subdirf2(fileres,"stde"),subdirf2(fileres,"stde")); + fprintf(fichtm,"\ + - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences (i.e population based) (estepm=%d months): %s
        \n", estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v")); fprintf(fichtm,"\ - - Health expectancies with their variances (no covariance): %s
        \n", + - Total life expectancy and total health expectancies to be spent in each health state e.j with their standard errors: %s
        \n", subdirf2(fileres,"t"),subdirf2(fileres,"t")); fprintf(fichtm,"\ - - Standard deviation of stable prevalences: %s
        \n",\ + - Standard deviation of period (stable) prevalences: %s
        \n",\ subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl")); /* if(popforecast==1) fprintf(fichtm,"\n */ @@ -3366,7 +3610,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); else fprintf(ficgp," \%%*lf (\%%*lf)"); } - fprintf(ficgp,"\" t\"Stable prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); + fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); else fprintf(ficgp," \%%*lf (\%%*lf)"); @@ -3414,7 +3658,8 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u for (k1=1; k1<= m ; k1 ++) { for (cpt=1; cpt<= nlstate ; cpt ++) { - k=2+nlstate*(2*cpt-2); + /* k=2+nlstate*(2*cpt-2); */ + k=2+(nlstate+1)*(cpt-1); fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1); fprintf(ficgp,"set ter png small\n\ set size 0.65,0.65\n\ @@ -3428,9 +3673,11 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u */ for (i=1; i< nlstate ; i ++) { - fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1); + fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+i,cpt,i+1); + /* fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/ } + fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+nlstate,cpt); } } @@ -3989,6 +4236,7 @@ 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;*/ @@ -4001,14 +4249,15 @@ double gompertz(double x[]) for (i=1;i<=imx ; i++) { - if (cens[i]==1 & wav[i]>1) + if (cens[i] == 1 && wav[i]>1) A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); - if (cens[i]==0 & wav[i]>1) + 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]);*/ } @@ -4069,6 +4318,7 @@ void printinggnuplotmort(char fileres[], + /***********************************************/ /**************** Main Program *****************/ /***********************************************/ @@ -4077,12 +4327,14 @@ 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; @@ -4099,6 +4351,7 @@ int main(int argc, char *argv[]) char line[MAXLINE], linepar[MAXLINE]; char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; char pathr[MAXLINE], pathimach[MAXLINE]; + char **bp, *tok, *val; /* pathtot */ int firstobs=1, lastobs=10; int sdeb, sfin; /* Status at beginning and end */ int c, h , cpt,l; @@ -4134,7 +4387,7 @@ int main(int argc, char *argv[]) char z[1]="c", occ; char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80]; - char strstart[80], *strt, strtend[80]; + char *strt, strtend[80]; char *stratrunc; int lstra; @@ -4176,7 +4429,17 @@ int main(int argc, char *argv[]) printf("\n%s\n%s",version,fullversion); if(argc <=1){ printf("\nEnter the parameter file name: "); - scanf("%s",pathtot); + fgets(pathr,FILENAMELENGTH,stdin); + i=strlen(pathr); + if(pathr[i-1]=='\n') + pathr[i-1]='\0'; + for (tok = pathr; tok != NULL; ){ + printf("Pathr |%s|\n",pathr); + while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0'); + printf("val= |%s| pathr=%s\n",val,pathr); + strcpy (pathtot, val); + if(pathr[0] == '\0') break; /* Dirty */ + } } else{ strcpy(pathtot,argv[1]); @@ -4194,7 +4457,9 @@ int main(int argc, char *argv[]) /* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */ split(pathtot,path,optionfile,optionfilext,optionfilefiname); printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); - chdir(path); + chdir(path); /* Can be a relative path */ + if(getcwd(pathcd,MAXLINE) > 0) /* So pathcd is the full path */ + printf("Current directory %s!\n",pathcd); strcpy(command,"mkdir "); strcat(command,optionfilefiname); if((outcmd=system(command)) != 0){ @@ -4360,7 +4625,6 @@ int main(int argc, char *argv[]) } fflush(ficlog); - p=param[1][1]; /* Reads comments: lines beginning with '#' */ @@ -4461,8 +4725,8 @@ int main(int argc, char *argv[]) /*-------- data file ----------*/ if((fic=fopen(datafile,"r"))==NULL) { - printf("Problem with datafile: %s\n", datafile);goto end; - fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end; + printf("Problem while opening datafile: %s\n", datafile);goto end; + fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);goto end; } n= lastobs; @@ -4484,43 +4748,125 @@ 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=0; line[j] != '\n';j++){ /* Untabifies line */ - if(line[j] == '\t') - line[j] = ' '; - } - for (j=maxwav;j>=1;j--){ - cutv(stra, strb,line,' '); s[j][i]=atoi(strb); - strcpy(line,stra); - 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; + dval=strtod(strb,&endptr); + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei); + exit(1); + } + weight[i]=dval; + 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 for individual %d, '%s'\n \ + Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ + for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ + For example, for multinomial values like 1, 2 and 3,\n \ + build V1=0 V2=0 for the reference value (1),\n \ + V1=1 V2=0 for (2) \n \ + and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ + output of IMaCh is often meaningless.\n \ + 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 */ + fclose(fic); /* printf("ii=%d", ij); scanf("%d",i);*/ imx=i-1; /* Number of individuals */ @@ -4618,8 +4964,7 @@ int main(int argc, char *argv[]) printf("cptcovprod=%d ", cptcovprod); fprintf(ficlog,"cptcovprod=%d ", cptcovprod); - scanf("%d ",i); - fclose(fic);*/ + scanf("%d ",i);*/ /* if(mle==1){*/ if (weightopt != 1) { /* Maximisation without weights*/ @@ -4810,11 +5155,13 @@ Title=%s
        Datafile=%s Firstpass=%d La \n\
        \

        • Parameter files

          \n\ + - Parameter file: %s.%s
          \n\ - Copy of the parameter file: o%s
          \n\ - Log file of the run: %s
          \n\ - Gnuplot file name: %s
          \n\ - Date and time at start: %s
        \n",\ fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\ + optionfilefiname,optionfilext,optionfilefiname,optionfilext,\ fileres,fileres,\ filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart); fflush(fichtm); @@ -4844,6 +5191,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); @@ -4853,9 +5201,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; } @@ -4864,12 +5212,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; } @@ -4878,29 +5230,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++) @@ -4918,48 +5270,48 @@ Interval (in months) between two waves: 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); + lsurv=vector(1,AGESUP); lpop=vector(1,AGESUP); tpop=vector(1,AGESUP); lsurv[agegomp]=100000; - - for (k=agegomp;k<=AGESUP;k++) { + + 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 */ @@ -5185,7 +5537,7 @@ lsurv=vector(1,AGESUP); /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);*/ /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ - replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */ + replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\ @@ -5207,18 +5559,18 @@ lsurv=vector(1,AGESUP); fclose(ficres); - /*--------------- Prevalence limit (stable prevalence) --------------*/ + /*--------------- Prevalence limit (period or stable prevalence) --------------*/ strcpy(filerespl,"pl"); strcat(filerespl,fileres); if((ficrespl=fopen(filerespl,"w"))==NULL) { - printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end; - fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end; + printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end; + fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end; } - printf("Computing 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"); + printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl); + fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl); + pstamp(ficrespl); + fprintf(ficrespl,"# Period (stable) prevalence \n"); fprintf(ficrespl,"#Age "); for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); fprintf(ficrespl,"\n"); @@ -5278,7 +5630,7 @@ lsurv=vector(1,AGESUP); hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ /* hstepm=1; aff par mois*/ - fprintf(ficrespij, "#Local time at start: %s", strstart); + pstamp(ficrespij); fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x "); for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ @@ -5348,8 +5700,8 @@ lsurv=vector(1,AGESUP); printf("Problem with total LE resultfile: %s\n", filerest);goto end; fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end; } - printf("Computing Total LEs with variances: file '%s' \n", filerest); - fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest); + printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); + fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); strcpy(filerese,"e"); @@ -5361,6 +5713,24 @@ lsurv=vector(1,AGESUP); printf("Computing Health Expectancies: result on file '%s' \n", filerese); fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese); + strcpy(fileresstde,"stde"); + strcat(fileresstde,fileres); + if((ficresstdeij=fopen(fileresstde,"w"))==NULL) { + printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0); + fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0); + } + printf("Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); + fprintf(ficlog,"Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde); + + strcpy(filerescve,"cve"); + strcat(filerescve,fileres); + if((ficrescveij=fopen(filerescve,"w"))==NULL) { + printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0); + fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0); + } + printf("Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); + fprintf(ficlog,"Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve); + strcpy(fileresv,"v"); strcat(fileresv,fileres); if((ficresvij=fopen(fileresv,"w"))==NULL) { @@ -5393,9 +5763,16 @@ lsurv=vector(1,AGESUP); fprintf(ficrest,"******\n"); fprintf(ficreseij,"\n#****** "); - for(j=1;j<=cptcoveff;j++) + fprintf(ficresstdeij,"\n#****** "); + fprintf(ficrescveij,"\n#****** "); + for(j=1;j<=cptcoveff;j++) { fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + } fprintf(ficreseij,"******\n"); + fprintf(ficresstdeij,"******\n"); + fprintf(ficrescveij,"******\n"); fprintf(ficresvij,"\n#****** "); for(j=1;j<=cptcoveff;j++) @@ -5404,7 +5781,8 @@ lsurv=vector(1,AGESUP); 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); + evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart); + cvevsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart); vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); oldm=oldms;savm=savms; @@ -5413,8 +5791,8 @@ lsurv=vector(1,AGESUP); 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) "); + pstamp(ficrest); + fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# Age ( e.. (std) "); for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); fprintf(ficrest,"\n"); @@ -5462,19 +5840,21 @@ lsurv=vector(1,AGESUP); free_ivector(cod,1,n); free_ivector(tab,1,NCOVMAX); fclose(ficreseij); + fclose(ficresstdeij); + fclose(ficrescveij); fclose(ficresvij); fclose(ficrest); fclose(ficpar); - /*------- Variance of stable prevalence------*/ + /*------- Variance of period (stable) prevalence------*/ strcpy(fileresvpl,"vpl"); strcat(fileresvpl,fileres); if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { - printf("Problem with variance of stable prevalence resultfile: %s\n", fileresvpl); + printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl); exit(0); } - printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl); + printf("Computing Variance-covariance of period (stable) prevalence: file '%s' \n", fileresvpl); for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ @@ -5498,11 +5878,11 @@ lsurv=vector(1,AGESUP); free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); } /* mle==-3 arrives here for freeing */ + free_matrix(prlim,1,nlstate,1,nlstate); free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); - free_matrix(covar,0,NCOVMAX,1,n); free_matrix(matcov,1,npar,1,npar); /*free_vector(delti,1,npar);*/ @@ -5517,7 +5897,8 @@ lsurv=vector(1,AGESUP); free_ivector(Tage,1,15); free_ivector(Tcode,1,100); - + free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX); + free_imatrix(codtab,1,100,1,10); fflush(fichtm); fflush(ficgp); @@ -5551,33 +5932,35 @@ lsurv=vector(1,AGESUP); fclose(ficlog); /*------ End -----------*/ - chdir(path); -#ifndef UNIX - /* strcpy(plotcmd,"\""); */ -#endif - strcpy(plotcmd,pathimach); + + printf("Before Current directory %s!\n",pathcd); + if(chdir(pathcd) != 0) + printf("Can't move to directory %s!\n",path); + if(getcwd(pathcd,MAXLINE) > 0) + printf("Current directory %s!\n",pathcd); /*strcat(plotcmd,CHARSEPARATOR);*/ - strcat(plotcmd,GNUPLOTPROGRAM); + sprintf(plotcmd,"gnuplot"); #ifndef UNIX - strcat(plotcmd,".exe"); - /* strcat(plotcmd,"\"");*/ + sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach); #endif - if(stat(plotcmd,&info)){ + if(!stat(plotcmd,&info)){ printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout); - } - -#ifndef UNIX - strcpy(plotcmd,"\""); -#endif - strcat(plotcmd,pathimach); - strcat(plotcmd,GNUPLOTPROGRAM); -#ifndef UNIX - strcat(plotcmd,".exe"); - strcat(plotcmd,"\""); + 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 - strcat(plotcmd," "); - strcat(plotcmd,optionfilegnuplot); - printf("Starting graphs with: %s",plotcmd);fflush(stdout); + }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("\n Problem with gnuplot\n");