--- imach/src/imach.c 2015/09/09 16:53:55 1.200 +++ imach/src/imach.c 2015/12/16 06:57:54 1.214 @@ -1,6 +1,64 @@ -/* $Id: imach.c,v 1.200 2015/09/09 16:53:55 brouard Exp $ +/* $Id: imach.c,v 1.214 2015/12/16 06:57:54 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.214 2015/12/16 06:57:54 brouard + Summary: temporary not working + + Revision 1.213 2015/12/11 18:22:17 brouard + Summary: 0.98r4 + + Revision 1.212 2015/11/21 12:47:24 brouard + Summary: minor typo + + Revision 1.211 2015/11/21 12:41:11 brouard + Summary: 0.98r3 with some graph of projected cross-sectional + + Author: Nicolas Brouard + + Revision 1.210 2015/11/18 17:41:20 brouard + Summary: Start working on projected prevalences + + Revision 1.209 2015/11/17 22:12:03 brouard + Summary: Adding ftolpl parameter + Author: N Brouard + + We had difficulties to get smoothed confidence intervals. It was due + to the period prevalence which wasn't computed accurately. The inner + parameter ftolpl is now an outer parameter of the .imach parameter + file after estepm. If ftolpl is small 1.e-4 and estepm too, + computation are long. + + Revision 1.208 2015/11/17 14:31:57 brouard + Summary: temporary + + Revision 1.207 2015/10/27 17:36:57 brouard + *** empty log message *** + + Revision 1.206 2015/10/24 07:14:11 brouard + *** empty log message *** + + Revision 1.205 2015/10/23 15:50:53 brouard + Summary: 0.98r3 some clarification for graphs on likelihood contributions + + Revision 1.204 2015/10/01 16:20:26 brouard + Summary: Some new graphs of contribution to likelihood + + Revision 1.203 2015/09/30 17:45:14 brouard + Summary: looking at better estimation of the hessian + + Also a better criteria for convergence to the period prevalence And + therefore adding the number of years needed to converge. (The + prevalence in any alive state shold sum to one + + Revision 1.202 2015/09/22 19:45:16 brouard + Summary: Adding some overall graph on contribution to likelihood. Might change + + Revision 1.201 2015/09/15 17:34:58 brouard + Summary: 0.98r0 + + - Some new graphs like suvival functions + - Some bugs fixed like model=1+age+V2. + Revision 1.200 2015/09/09 16:53:55 brouard Summary: Big bug thanks to Flavia @@ -638,6 +696,10 @@ /* #define DEBUG */ /* #define DEBUGBRENT */ +/* #define DEBUGLINMIN */ +/* #define DEBUGHESS */ +#define DEBUGHESSIJ +/* #define LINMINORIGINAL /\* Don't use loop on scale in linmin (accepting nan)*\/ */ #define POWELL /* Instead of NLOPT */ #define POWELLF1F3 /* Skip test */ /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */ @@ -709,6 +771,8 @@ typedef struct { #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 +/*#define decodtabm(h,k,cptcoveff)= (h <= (1<> (k-1)) & 1) +1 : -1)*/ +#define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 #define MAXN 20000 #define YEARM 12. /**< Number of months per year */ #define AGESUP 130 @@ -725,12 +789,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.200 2015/09/09 16:53:55 brouard Exp $ */ +/* $Id: imach.c,v 1.214 2015/12/16 06:57:54 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; -char copyright[]="September 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; -char fullversion[]="$Revision: 1.200 $ $Date: 2015/09/09 16:53:55 $"; +char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; +char fullversion[]="$Revision: 1.214 $ $Date: 2015/12/16 06:57:54 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -766,7 +830,7 @@ double **matprod2(); /* test */ double **oldm, **newm, **savm; /* Working pointers to matrices */ double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ /*FILE *fic ; */ /* Used in readdata only */ -FILE *ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; +FILE *ficpar, *ficparo,*ficres, *ficresp, *ficresphtm, *ficresphtmfr, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; FILE *ficlog, *ficrespow; int globpr=0; /* Global variable for printing or not */ double fretone; /* Only one call to likelihood */ @@ -796,7 +860,7 @@ char command[FILENAMELENGTH]; int outcmd=0; char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH]; - +char fileresu[FILENAMELENGTH]; /* fileres without r in front */ char filelog[FILENAMELENGTH]; /* Log file */ char filerest[FILENAMELENGTH]; char fileregp[FILENAMELENGTH]; @@ -920,7 +984,7 @@ static int split( char *path, char *dirc } /* got dirc from getcwd*/ printf(" DIRC = %s \n",dirc); - } else { /* strip direcotry from path */ + } else { /* strip directory from path */ ss++; /* after this, the filename */ l2 = strlen( ss ); /* length of filename */ if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH ); @@ -1328,7 +1392,30 @@ char *subdirf3(char fileres[], char *pre strcat(tmpout,fileres); return tmpout; } + +/*************** function subdirfext ***********/ +char *subdirfext(char fileres[], char *preop, char *postop) +{ + + strcpy(tmpout,preop); + strcat(tmpout,fileres); + strcat(tmpout,postop); + return tmpout; +} +/*************** function subdirfext3 ***********/ +char *subdirfext3(char fileres[], char *preop, char *postop) +{ + + /* Caution optionfilefiname is hidden */ + strcpy(tmpout,optionfilefiname); + strcat(tmpout,"/"); + strcat(tmpout,preop); + strcat(tmpout,fileres); + strcat(tmpout,postop); + return tmpout; +} + char *asc_diff_time(long time_sec, char ascdiff[]) { long sec_left, days, hours, minutes; @@ -1592,23 +1679,29 @@ void linmin(double p[], double xi[], int double xx,xmin,bx,ax; double fx,fb,fa; - double scale=10., axs, xxs, xxss; /* Scale added for infinity */ - +#ifdef LINMINORIGINAL +#else + double scale=10., axs, xxs; /* Scale added for infinity */ +#endif + ncom=n; pcom=vector(1,n); xicom=vector(1,n); nrfunc=func; for (j=1;j<=n;j++) { pcom[j]=p[j]; - xicom[j]=xi[j]; + xicom[j]=xi[j]; /* Former scale xi[j] of currrent direction i */ } - /* axs=0.0; */ - /* xxss=1; /\* 1 and using scale *\/ */ - xxs=1; - /* do{ */ - ax=0.; +#ifdef LINMINORIGINAL + xx=1.; +#else + axs=0.0; + xxs=1.; + do{ xx= xxs; +#endif + ax=0.; mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Outputs: xtx[j]=pcom[j]+(*xx)*xicom[j]; fx=f(xtx[j]) */ /* brackets with inputs ax=0 and xx=1, but points, pcom=p, and directions values, xicom=xi, are sent via f1dim(x) */ /* xt[x,j]=pcom[j]+x*xicom[j] f(ax) = f(xt(a,j=1,n)) = f(p(j) + 0 * xi(j)) and f(xx) = f(xt(x, j=1,n)) = f(p(j) + 1 * xi(j)) */ @@ -1616,14 +1709,22 @@ void linmin(double p[], double xi[], int /* Given input ax=axs and xx=xxs, xx might be too far from ax to get a finite f(xx) */ /* Searches on line, outputs (ax, xx, bx) such that fx < min(fa and fb) */ /* Find a bracket a,x,b in direction n=xi ie xicom, order may change. Scale is [0:xxs*xi[j]] et non plus [0:xi[j]]*/ - /* if (fx != fx){ */ - /* xxs=xxs/scale; /\* Trying a smaller xx, closer to initial ax=0 *\/ */ - /* printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx); */ - /* } */ - /* }while(fx != fx); */ - +#ifdef LINMINORIGINAL +#else + if (fx != fx){ + xxs=xxs/scale; /* Trying a smaller xx, closer to initial ax=0 */ + printf("|"); + fprintf(ficlog,"|"); +#ifdef DEBUGLINMIN + printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx); +#endif + } + }while(fx != fx); +#endif + #ifdef DEBUGLINMIN printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb); + fprintf(ficlog,"\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb); #endif *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Giving a bracketting triplet (ax, xx, bx), find a minimum, xmin, according to f1dim, *fret(xmin),*/ /* fa = f(p[j] + ax * xi[j]), fx = f(p[j] + xx * xi[j]), fb = f(p[j] + bx * xi[j]) */ @@ -1636,22 +1737,37 @@ void linmin(double p[], double xi[], int #endif #ifdef DEBUGLINMIN printf("linmin end "); + fprintf(ficlog,"linmin end "); #endif for (j=1;j<=n;j++) { - /* printf(" before xi[%d]=%12.8f", j,xi[j]); */ - xi[j] *= xmin; /* xi rescaled by xmin: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */ - /* if(xxs <1.0) */ - /* printf(" after xi[%d]=%12.8f, xmin=%12.8f, ax=%12.8f, xx=%12.8f, bx=%12.8f, xxs=%12.8f", j,xi[j], xmin, ax, xx, bx,xxs ); */ +#ifdef LINMINORIGINAL + xi[j] *= xmin; +#else +#ifdef DEBUGLINMIN + if(xxs <1.0) + printf(" before xi[%d]=%12.8f", j,xi[j]); +#endif + xi[j] *= xmin*xxs; /* xi rescaled by xmin and number of loops: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */ +#ifdef DEBUGLINMIN + if(xxs <1.0) + printf(" after xi[%d]=%12.8f, xmin=%12.8f, ax=%12.8f, xx=%12.8f, bx=%12.8f, xxs=%12.8f", j,xi[j], xmin, ax, xx, bx,xxs ); +#endif +#endif p[j] += xi[j]; /* Parameters values are updated accordingly */ } - /* printf("\n"); */ #ifdef DEBUGLINMIN + printf("\n"); printf("Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p)); + fprintf(ficlog,"Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p)); for (j=1;j<=n;j++) { - printf(" xi[%d]= %12.7f p[%d]= %12.7f",j,xi[j],j,p[j]); - if(j % ncovmodel == 0) + printf(" xi[%d]= %14.10f p[%d]= %12.7f",j,xi[j],j,p[j]); + fprintf(ficlog," xi[%d]= %14.10f p[%d]= %12.7f",j,xi[j],j,p[j]); + if(j % ncovmodel == 0){ printf("\n"); + fprintf(ficlog,"\n"); + } } +#else #endif free_vector(xicom,1,n); free_vector(pcom,1,n); @@ -1685,7 +1801,7 @@ void powell(double p[], double **xi, int xits=vector(1,n); *fret=(*func)(p); for (j=1;j<=n;j++) pt[j]=p[j]; - rcurr_time = time(NULL); + rcurr_time = time(NULL); for (*iter=1;;++(*iter)) { fp=(*fret); /* From former iteration or initial value */ ibig=0; @@ -1729,10 +1845,10 @@ void powell(double p[], double **xi, int for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */ fptt=(*fret); #ifdef DEBUG - printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret); - fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret); + printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret); + fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret); #endif - printf("%d",i);fflush(stdout); /* print direction (parameter) i */ + printf("%d",i);fflush(stdout); /* print direction (parameter) i */ fprintf(ficlog,"%d",i);fflush(ficlog); linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input.*/ /* Outputs are fret(new point p) p is updated and xit rescaled */ @@ -1830,7 +1946,7 @@ void powell(double p[], double **xi, int t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); /* Intel compiler doesn't work on one line; bug reported */ t= t- del*SQR(fp-fptt); #endif - directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If del was big enough we change it for a new direction */ + directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If delta was big enough we change it for a new direction */ #ifdef DEBUG printf("t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest); fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest); @@ -1845,9 +1961,9 @@ void powell(double p[], double **xi, int if (t < 0.0) { /* Then we use it for new direction */ #else if (directest*t < 0.0) { /* Contradiction between both tests */ - printf("directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del); + printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del); printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); - fprintf(ficlog,"directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del); + fprintf(ficlog,"directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del); fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); } if (directest < 0.0) { /* Then we use it for new direction */ @@ -1855,17 +1971,23 @@ void powell(double p[], double **xi, int #ifdef DEBUGLINMIN printf("Before linmin in direction P%d-P0\n",n); for (j=1;j<=n;j++) { - printf("Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); - if(j % ncovmodel == 0) + printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + if(j % ncovmodel == 0){ printf("\n"); + fprintf(ficlog,"\n"); + } } #endif linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/ #ifdef DEBUGLINMIN for (j=1;j<=n;j++) { printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); - if(j % ncovmodel == 0) + fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); + if(j % ncovmodel == 0){ printf("\n"); + fprintf(ficlog,"\n"); + } } #endif for (j=1;j<=n;j++) { @@ -1895,18 +2017,40 @@ void powell(double p[], double **xi, int /**** Prevalence limit (stable or period prevalence) ****************/ -double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij) +double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij) { /* Computes the prevalence limit in each live state at age x by left multiplying the unit - matrix by transitions matrix until convergence is reached */ - + matrix by transitions matrix until convergence is reached with precision ftolpl */ + /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */ + /* Wx is row vector: population in state 1, population in state 2, population dead */ + /* or prevalence in state 1, prevalence in state 2, 0 */ + /* newm is the matrix after multiplications, its rows are identical at a factor */ + /* Initial matrix pimij */ + /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */ + /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */ + /* 0, 0 , 1} */ + /* + * and after some iteration: */ + /* {0.45504275246439968, 0.42731458730878791, 0.11764266022681241, */ + /* 0.45201005341706885, 0.42865420071559901, 0.11933574586733192, */ + /* 0, 0 , 1} */ + /* And prevalence by suppressing the deaths are close to identical rows in prlim: */ + /* {0.51571254859325999, 0.4842874514067399, */ + /* 0.51326036147820708, 0.48673963852179264} */ + /* If we start from prlim again, prlim tends to a constant matrix */ + int i, ii,j,k; - double min, max, maxmin, maxmax,sumnew=0.; + double *min, *max, *meandiff, maxmax,sumnew=0.; /* double **matprod2(); */ /* test */ double **out, cov[NCOVMAX+1], **pmij(); double **newm; - double agefin, delaymax=50 ; /* Max number of years to converge */ + double agefin, delaymax=200. ; /* 100 Max number of years to converge */ + int ncvloop=0; + min=vector(1,nlstate); + max=vector(1,nlstate); + meandiff=vector(1,nlstate); + for (ii=1;ii<=nlstate+ndeath;ii++) for (j=1;j<=nlstate+ndeath;j++){ oldm[ii][j]=(ii==j ? 1.0 : 0.0); @@ -1915,7 +2059,9 @@ double **prevalim(double **prlim, int nl cov[1]=1.; /* Even if hstepm = 1, at least one multiplication by the unit matrix */ + /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */ for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){ + ncvloop++; newm=savm; /* Covariates have to be included here again */ cov[2]=agefin; @@ -1942,25 +2088,45 @@ double **prevalim(double **prlim, int nl savm=oldm; oldm=newm; - maxmax=0.; - for(j=1;j<=nlstate;j++){ - min=1.; - max=0.; - for(i=1; i<=nlstate; i++) { - sumnew=0; - for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; + + for(j=1; j<=nlstate; j++){ + max[j]=0.; + min[j]=1.; + } + for(i=1;i<=nlstate;i++){ + sumnew=0; + for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; + for(j=1; j<=nlstate; j++){ prlim[i][j]= newm[i][j]/(1-sumnew); - /*printf(" prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d \n", i, j, i, j, prlim[i][j],(int)agefin);*/ - max=FMAX(max,prlim[i][j]); - min=FMIN(min,prlim[i][j]); + max[j]=FMAX(max[j],prlim[i][j]); + min[j]=FMIN(min[j],prlim[i][j]); } - maxmin=max-min; - maxmax=FMAX(maxmax,maxmin); + } + + maxmax=0.; + for(j=1; j<=nlstate; j++){ + meandiff[j]=(max[j]-min[j])/(max[j]+min[j])*2.; /* mean difference for each column */ + maxmax=FMAX(maxmax,meandiff[j]); + /* printf(" age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, j, meandiff[j],(int)agefin, j, max[j], j, min[j],maxmax); */ } /* j loop */ + *ncvyear= (int)age- (int)agefin; + /* printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */ if(maxmax < ftolpl){ + /* printf("maxmax=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */ + free_vector(min,1,nlstate); + free_vector(max,1,nlstate); + free_vector(meandiff,1,nlstate); return prlim; } } /* age loop */ + /* After some age loop it doesn't converge */ + printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ +Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); + /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ + free_vector(min,1,nlstate); + free_vector(max,1,nlstate); + free_vector(meandiff,1,nlstate); + return prlim; /* should not reach here */ } @@ -2086,6 +2252,7 @@ double ***hpxij(double ***po, int nhstep double **out, cov[NCOVMAX+1]; double **newm; double agexact; + double agebegin, ageend; /* Hstepm could be zero and should return the unit matrix */ for (i=1;i<=nlstate+ndeath;i++) @@ -2099,7 +2266,7 @@ double ***hpxij(double ***po, int nhstep newm=savm; /* Covariates have to be included here again */ cov[1]=1.; - agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; + agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */ cov[2]=agexact; if(nagesqr==1) cov[3]= agexact*agexact; @@ -2470,6 +2637,7 @@ double funcone( double *x) int s1, s2; double bbh, survp; double agexact; + double agebegin, ageend; /*extern weight */ /* We are differentiating ll according to initial status */ /* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ @@ -2488,7 +2656,12 @@ double funcone( double *x) oldm[ii][j]=(ii==j ? 1.0 : 0.0); savm[ii][j]=(ii==j ? 1.0 : 0.0); } - for(d=0; dFile of contributions to the likelihood: %s
\n",subdirf(fileresilk),subdirf(fileresilk)); - fflush(fichtm); - } + if (mle ==0) + fprintf(fichtm,"\n
File of contributions to the likelihood computed with initial parameters and mle = %d.",mle); + else if(mle >=1) + fprintf(fichtm,"\n
File of contributions to the likelihood computed with optimized parameters mle = %d.",mle); + fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: %s
\n",subdirf(fileresilk),subdirf(fileresilk)); + + + for (k=1; k<= nlstate ; k++) { + fprintf(fichtm,"
- Probability p%dj by origin %d and destination j. Dot's sizes are related to corresponding weight: %s-p%dj.png
\ +",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); + } + fprintf(fichtm,"
- The function drawn is -2Log(L) in Log scale: by state of origin %s-ori.png
\ +",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); + fprintf(fichtm,"
- and by state of destination %s-dest.png
\ +",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); + fflush(fichtm); + } return; } @@ -2620,7 +2807,7 @@ void mlikeli(FILE *ficres,double p[], in for (j=1;j<=npar;j++) xi[i][j]=(i==j ? 1.0 : 0.0); printf("Powell\n"); fprintf(ficlog,"Powell\n"); - strcpy(filerespow,"pow"); + strcpy(filerespow,"POW_"); strcat(filerespow,fileres); if((ficrespow=fopen(filerespow,"w"))==NULL) { printf("Problem with resultfile: %s\n", filerespow); @@ -2663,32 +2850,32 @@ void mlikeli(FILE *ficres,double p[], in #endif free_matrix(xi,1,npar,1,npar); fclose(ficrespow); - printf("#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); - fprintf(ficlog,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); + printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); + fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); fprintf(ficres,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); } /**** Computes Hessian and covariance matrix ***/ -void hesscov(double **matcov, double p[], int npar, double delti[], double ftolhess, double (*func)(double [])) +void hesscov(double **matcov, double **hess, double p[], int npar, double delti[], double ftolhess, double (*func)(double [])) { double **a,**y,*x,pd; - double **hess; + /* double **hess; */ int i, j; int *indx; double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar); - double hessij(double p[], double delti[], int i, int j,double (*func)(double []),int npar); + double hessij(double p[], double **hess, double delti[], int i, int j,double (*func)(double []),int npar); void lubksb(double **a, int npar, int *indx, double b[]) ; void ludcmp(double **a, int npar, int *indx, double *d) ; double gompertz(double p[]); - hess=matrix(1,npar,1,npar); + /* hess=matrix(1,npar,1,npar); */ printf("\nCalculation of the hessian matrix. Wait...\n"); fprintf(ficlog,"\nCalculation of the hessian matrix. Wait...\n"); for (i=1;i<=npar;i++){ - printf("%d",i);fflush(stdout); - fprintf(ficlog,"%d",i);fflush(ficlog); + printf("%d-",i);fflush(stdout); + fprintf(ficlog,"%d-",i);fflush(ficlog); hess[i][i]=hessii(p,ftolhess,i,delti,func,npar); @@ -2699,9 +2886,9 @@ void hesscov(double **matcov, double p[] for (i=1;i<=npar;i++) { for (j=1;j<=npar;j++) { if (j>i) { - printf(".%d%d",i,j);fflush(stdout); - fprintf(ficlog,".%d%d",i,j);fflush(ficlog); - hess[i][j]=hessij(p,delti,i,j,func,npar); + printf(".%d-%d",i,j);fflush(stdout); + fprintf(ficlog,".%d-%d",i,j);fflush(ficlog); + hess[i][j]=hessij(p,hess, delti,i,j,func,npar); hess[j][i]=hess[i][j]; /*printf(" %lf ",hess[i][j]);*/ @@ -2735,53 +2922,78 @@ void hesscov(double **matcov, double p[] fprintf(ficlog,"\n#Hessian matrix#\n"); for (i=1;i<=npar;i++) { for (j=1;j<=npar;j++) { - printf("%.3e ",hess[i][j]); - fprintf(ficlog,"%.3e ",hess[i][j]); + printf("%.6e ",hess[i][j]); + fprintf(ficlog,"%.6e ",hess[i][j]); } printf("\n"); fprintf(ficlog,"\n"); } + /* printf("\n#Covariance matrix#\n"); */ + /* fprintf(ficlog,"\n#Covariance matrix#\n"); */ + /* for (i=1;i<=npar;i++) { */ + /* for (j=1;j<=npar;j++) { */ + /* printf("%.6e ",matcov[i][j]); */ + /* fprintf(ficlog,"%.6e ",matcov[i][j]); */ + /* } */ + /* printf("\n"); */ + /* fprintf(ficlog,"\n"); */ + /* } */ + /* Recompute Inverse */ - for (i=1;i<=npar;i++) - for (j=1;j<=npar;j++) a[i][j]=matcov[i][j]; - ludcmp(a,npar,indx,&pd); + /* for (i=1;i<=npar;i++) */ + /* for (j=1;j<=npar;j++) a[i][j]=matcov[i][j]; */ + /* ludcmp(a,npar,indx,&pd); */ + + /* printf("\n#Hessian matrix recomputed#\n"); */ + + /* for (j=1;j<=npar;j++) { */ + /* for (i=1;i<=npar;i++) x[i]=0; */ + /* x[j]=1; */ + /* lubksb(a,npar,indx,x); */ + /* for (i=1;i<=npar;i++){ */ + /* y[i][j]=x[i]; */ + /* printf("%.3e ",y[i][j]); */ + /* fprintf(ficlog,"%.3e ",y[i][j]); */ + /* } */ + /* printf("\n"); */ + /* fprintf(ficlog,"\n"); */ + /* } */ - /* printf("\n#Hessian matrix recomputed#\n"); + /* Verifying the inverse matrix */ +#ifdef DEBUGHESS + y=matprod2(y,hess,1,npar,1,npar,1,npar,matcov); + + printf("\n#Verification: multiplying the matrix of covariance by the Hessian matrix, should be unity:#\n"); + fprintf(ficlog,"\n#Verification: multiplying the matrix of covariance by the Hessian matrix. Should be unity:#\n"); for (j=1;j<=npar;j++) { - for (i=1;i<=npar;i++) x[i]=0; - x[j]=1; - lubksb(a,npar,indx,x); for (i=1;i<=npar;i++){ - y[i][j]=x[i]; - printf("%.3e ",y[i][j]); - fprintf(ficlog,"%.3e ",y[i][j]); + printf("%.2f ",y[i][j]); + fprintf(ficlog,"%.2f ",y[i][j]); } printf("\n"); fprintf(ficlog,"\n"); } - */ +#endif free_matrix(a,1,npar,1,npar); free_matrix(y,1,npar,1,npar); free_vector(x,1,npar); free_ivector(indx,1,npar); - free_matrix(hess,1,npar,1,npar); + /* free_matrix(hess,1,npar,1,npar); */ } /*************** hessian matrix ****************/ double hessii(double x[], double delta, int theta, double delti[], double (*func)(double []), int npar) -{ +{ /* Around values of x, computes the function func and returns the scales delti and hessian */ int i; int l=1, lmax=20; - double k1,k2; + double k1,k2, res, fx; double p2[MAXPARM+1]; /* identical to x */ - double res; double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; - double fx; int k=0,kmax=10; double l1; @@ -2797,9 +3009,9 @@ double hessii(double x[], double delta, p2[theta]=x[theta]-delt; k2=func(p2)-fx; /*res= (k1-2.0*fx+k2)/delt/delt; */ - res= (k1+k2)/delt/delt/2.; /* Divided by because L and not 2*L */ + res= (k1+k2)/delt/delt/2.; /* Divided by 2 because L and not 2*L */ -#ifdef DEBUGHESS +#ifdef DEBUGHESSII printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); fprintf(ficlog,"%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); #endif @@ -2813,48 +3025,125 @@ double hessii(double x[], double delta, else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){ delts=delt; } - } + } /* End loop k */ } delti[theta]=delts; return res; } -double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) +double hessij( double x[], double **hess, double delti[], int thetai,int thetaj,double (*func)(double []),int npar) { int i; int l=1, lmax=20; double k1,k2,k3,k4,res,fx; double p2[MAXPARM+1]; - int k; + int k, kmax=1; + double v1, v2, cv12, lc1, lc2; + int firstime=0; + fx=func(x); - for (k=1; k<=2; k++) { + for (k=1; k<=kmax; k=k+10) { for (i=1;i<=npar;i++) p2[i]=x[i]; - p2[thetai]=x[thetai]+delti[thetai]/k; - p2[thetaj]=x[thetaj]+delti[thetaj]/k; + p2[thetai]=x[thetai]+delti[thetai]*k; + p2[thetaj]=x[thetaj]+delti[thetaj]*k; k1=func(p2)-fx; - p2[thetai]=x[thetai]+delti[thetai]/k; - p2[thetaj]=x[thetaj]-delti[thetaj]/k; + p2[thetai]=x[thetai]+delti[thetai]*k; + p2[thetaj]=x[thetaj]-delti[thetaj]*k; k2=func(p2)-fx; - p2[thetai]=x[thetai]-delti[thetai]/k; - p2[thetaj]=x[thetaj]+delti[thetaj]/k; + p2[thetai]=x[thetai]-delti[thetai]*k; + p2[thetaj]=x[thetaj]+delti[thetaj]*k; k3=func(p2)-fx; - p2[thetai]=x[thetai]-delti[thetai]/k; - p2[thetaj]=x[thetaj]-delti[thetaj]/k; + p2[thetai]=x[thetai]-delti[thetai]*k; + p2[thetaj]=x[thetaj]-delti[thetaj]*k; k4=func(p2)-fx; - res=(k1-k2-k3+k4)/4.0/delti[thetai]*k/delti[thetaj]*k/2.; /* Because of L not 2*L */ -#ifdef DEBUG - printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); - fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + res=(k1-k2-k3+k4)/4.0/delti[thetai]/k/delti[thetaj]/k/2.; /* Because of L not 2*L */ + if(k1*k2*k3*k4 <0.){ + firstime=1; + kmax=kmax+10; + } + if(kmax >=10 || firstime ==1){ + printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol); + fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol); + printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + } +#ifdef DEBUGHESSIJ + v1=hess[thetai][thetai]; + v2=hess[thetaj][thetaj]; + cv12=res; + /* Computing eigen value of Hessian matrix */ + lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; + lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; + if ((lc2 <0) || (lc1 <0) ){ + printf("Warning: sub Hessian matrix '%d%d' does not have positive eigen values \n",thetai,thetaj); + fprintf(ficlog, "Warning: sub Hessian matrix '%d%d' does not have positive eigen values \n",thetai,thetaj); + printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); + } #endif } return res; } + /* Not done yet: Was supposed to fix if not exactly at the maximum */ +/* double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) */ +/* { */ +/* int i; */ +/* int l=1, lmax=20; */ +/* double k1,k2,k3,k4,res,fx; */ +/* double p2[MAXPARM+1]; */ +/* double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; */ +/* int k=0,kmax=10; */ +/* double l1; */ + +/* fx=func(x); */ +/* for(l=0 ; l <=lmax; l++){ /\* Enlarging the zone around the Maximum *\/ */ +/* l1=pow(10,l); */ +/* delts=delt; */ +/* for(k=1 ; k khi/nkhif) || (k2 >khi/nkhif) || (k4 >khi/nkhif) || (k4 >khi/nkhif)){ /\* Keeps lastvalue before 3.84/2 KHI2 5% 1d.f. *\/ */ +/* k=kmax; l=lmax*10; */ +/* } */ +/* else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){ */ +/* delts=delt; */ +/* } */ +/* } /\* End loop k *\/ */ +/* } */ +/* delti[theta]=delts; */ +/* return res; */ +/* } */ + + /************** Inverse of matrix **************/ void ludcmp(double **a, int n, int *indx, double *d) { @@ -2935,25 +3224,61 @@ void pstamp(FILE *fichier) } /************ 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[]) +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[],\ + int firstpass, int lastpass, int stepm, int weightopt, char model[]) { /* Some frequencies */ int i, m, jk, j1, bool, z1,j; + int mi; /* Effective wave */ int first; double ***freq; /* Frequencies */ double *pp, **prop; double pos,posprop, k2, dateintsum=0,k2cpt=0; - char fileresp[FILENAMELENGTH]; - + char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH]; + double agebegin, ageend; + pp=vector(1,nlstate); prop=matrix(1,nlstate,iagemin,iagemax+3); - strcpy(fileresp,"p"); - strcat(fileresp,fileres); + strcpy(fileresp,"P_"); + strcat(fileresp,fileresu); + /*strcat(fileresphtm,fileresu);*/ if((ficresp=fopen(fileresp,"w"))==NULL) { printf("Problem with prevalence resultfile: %s\n", fileresp); fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); exit(0); } + + strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm")); + if((ficresphtm=fopen(fileresphtm,"w"))==NULL) { + printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); + fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); + fflush(ficlog); + exit(70); + } + else{ + fprintf(ficresphtm,"\nIMaCh PHTM_ %s\n %s
%s \ +
\n\ +Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\ + fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); + } + fprintf(ficresphtm,"Current page is file %s
\n\n

Frequencies and prevalence by age at begin of transition

\n",fileresphtm, fileresphtm); + + strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); + if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { + printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); + fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); + fflush(ficlog); + exit(70); + } + else{ + fprintf(ficresphtmfr,"\nIMaCh PHTM_Frequency table %s\n %s
%s \ +
\n\ +Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\ + fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); + } + fprintf(ficresphtmfr,"Current page is file %s
\n\n

Frequencies of all effective transitions by age at begin of transition

Unknown status is -1
\n",fileresphtmfr, fileresphtmfr); + freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3); j1=0; @@ -2962,10 +3287,7 @@ void freqsummary(char fileres[], int ia first=1; - /* for(k1=1; k1<=j ; k1++){ */ /* Loop on covariates */ - /* for(i1=1; i1<=ncodemax[k1];i1++){ */ /* Now it is 2 */ - /* j1++; */ - for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ + for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ /* Loop on covariates combination */ /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); scanf("%d", i);*/ for (i=-5; i<=nlstate+ndeath; i++) @@ -2979,9 +3301,9 @@ void freqsummary(char fileres[], int ia dateintsum=0; k2cpt=0; - for (i=1; i<=imx; i++) { + for (i=1; i<=imx; i++) { /* For each individual i */ bool=1; - if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ + if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ for (z1=1; z1<=cptcoveff; z1++) if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ /* Tests if the value of each of the covariates of i is equal to filter j1 */ @@ -2991,51 +3313,98 @@ void freqsummary(char fileres[], int ia j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ } - } - + } /* cptcovn > 0 */ + if (bool==1){ - for(m=firstpass; m<=lastpass; m++){ - k2=anint[m][i]+(mint[m][i]/12.); - /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ - if(agev[m][i]==0) agev[m][i]=iagemax+1; - if(agev[m][i]==1) agev[m][i]=iagemax+2; - if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i]; + /* for(m=firstpass; m<=lastpass; m++){ */ + for(mi=1; mi=firstpass && m <=lastpass){ + k2=anint[m][i]+(mint[m][i]/12.); + /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ + if(agev[m][i]==0) agev[m][i]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */ + if(agev[m][i]==1) agev[m][i]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */ + if (s[m][i]>0 && s[m][i]<=nlstate) /* If status at wave m is known and a live state */ + prop[s[m][i]][(int)agev[m][i]] += weight[i]; /* At age of beginning of transition, where status is known */ if (m1) && (agev[m][i]< (iagemax+3))) { - dateintsum=dateintsum+k2; - k2cpt++; + /* if(s[m][i]==4 && s[m+1][i]==4) */ + /* printf(" num=%ld m=%d, i=%d s1=%d s2=%d agev at m=%d\n", num[i], m, i,s[m][i],s[m+1][i], (int)agev[m][i]); */ + if(s[m][i]==-1) + printf(" num=%ld m=%d, i=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[i], m, i,s[m][i],s[m+1][i], (int)agev[m][i],agebegin, ageend, (int)((agebegin+ageend)/2.)); + freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; /* At age of beginning of transition, where status is known */ + /* freq[s[m][i]][s[m+1][i]][(int)((agebegin+ageend)/2.)] += weight[i]; */ + freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */ } - /*}*/ - } - } - } /* end i */ + } + if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3)) && (anint[m][i]!=9999) && (mint[m][i]!=99)) { + dateintsum=dateintsum+k2; + k2cpt++; + /* printf("i=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",i, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ + } + /*}*/ + } /* end m */ + } /* end bool */ + } /* end i = 1 to imx */ /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ pstamp(ficresp); if (cptcovn>0) { fprintf(ficresp, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresp, "**********\n#"); + fprintf(ficresphtm, "\n

********** Variable "); + fprintf(ficresphtmfr, "\n

********** Variable "); + for (z1=1; z1<=cptcoveff; z1++){ + fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + } + fprintf(ficresp, "**********\n#"); + fprintf(ficresphtm, "**********

\n"); + fprintf(ficresphtmfr, "**********\n"); fprintf(ficlog, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficlog, "**********\n#"); + fprintf(ficlog, "**********\n"); } - for(i=1; i<=nlstate;i++) + fprintf(ficresphtm,""); + for(i=1; i<=nlstate;i++) { fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); + fprintf(ficresphtm, "",i,i); + } fprintf(ficresp, "\n"); + fprintf(ficresphtm, "\n"); + + /* Header of frequency table by age */ + fprintf(ficresphtmfr,"
AgePrev(%d)N(%d)N
"); + fprintf(ficresphtmfr," "); + for(jk=-1; jk <=nlstate+ndeath; jk++){ + for(m=-1; m <=nlstate+ndeath; m++){ + if(jk!=0 && m!=0) + fprintf(ficresphtmfr," ",jk,m); + } + } + fprintf(ficresphtmfr, "\n"); + /* For each age */ for(i=iagemin; i <= iagemax+3; i++){ - if(i==iagemax+3){ + fprintf(ficresphtm,""); + if(i==iagemax+1){ + fprintf(ficlog,"1"); + fprintf(ficresphtmfr," "); + }else if(i==iagemax+2){ + fprintf(ficlog,"0"); + fprintf(ficresphtmfr," "); + }else if(i==iagemax+3){ fprintf(ficlog,"Total"); + fprintf(ficresphtmfr," "); }else{ if(first==1){ first=0; printf("See log file for details...\n"); } + fprintf(ficresphtmfr," ",i); fprintf(ficlog,"Age %d", i); } for(jk=1; jk <=nlstate ; jk++){ @@ -3078,32 +3447,47 @@ void freqsummary(char fileres[], int ia if( i <= iagemax){ if(pos>=1.e-5){ fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); + fprintf(ficresphtm,"",i,prop[jk][i]/posprop, prop[jk][i],posprop); /*probs[i][jk][j1]= pp[jk]/pos;*/ /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/ } - else + else{ fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop); + fprintf(ficresphtm,"",i, prop[jk][i],posprop); + } } } - for(jk=-1; jk <=nlstate+ndeath; jk++) - for(m=-1; m <=nlstate+ndeath; m++) - if(freq[jk][m][i] !=0 ) { - if(first==1) - printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); + for(jk=-1; jk <=nlstate+ndeath; jk++){ + for(m=-1; m <=nlstate+ndeath; m++){ + if(freq[jk][m][i] !=0 ) { /* minimizing output */ + if(first==1){ + printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); + } fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]); } - if(i <= iagemax) + if(jk!=0 && m!=0) + fprintf(ficresphtmfr," ",freq[jk][m][i]); + } + } + fprintf(ficresphtmfr,"\n "); + if(i <= iagemax){ fprintf(ficresp,"\n"); + fprintf(ficresphtm,"\n"); + } if(first==1) printf("Others in log...\n"); fprintf(ficlog,"\n"); - } + } /* end loop i */ + fprintf(ficresphtm,"
Age%d%d
0
Unknown
Total
%d%d%.5f%.0f%.0f%dNaNq%.0f%.0f%.0f
\n"); + fprintf(ficresphtmfr,"\n"); /*}*/ - } + } /* end j1 */ dateintmean=dateintsum/k2cpt; fclose(ficresp); + fclose(ficresphtm); + fclose(ficresphtmfr); 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); @@ -3119,6 +3503,9 @@ void prevalence(double ***probs, double */ int i, m, jk, j1, bool, z1,j; + int mi; /* Effective wave */ + int iage; + double agebegin, ageend; double **prop; double posprop; @@ -3138,54 +3525,57 @@ void prevalence(double ***probs, double first=1; for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ - /*for(i1=1; i1<=ncodemax[k1];i1++){ - j1++;*/ - - for (i=1; i<=nlstate; i++) - for(m=iagemin; m <= iagemax+3; m++) - prop[i][m]=0.0; - - for (i=1; i<=imx; i++) { /* Each individual */ - bool=1; - if (cptcovn>0) { - for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) - bool=0; - } - if (bool==1) { - for(m=firstpass; m<=lastpass; m++){/* Other selection (we can limit to certain interviews*/ + for (i=1; i<=nlstate; i++) + for(iage=iagemin; iage <= iagemax+3; iage++) + prop[i][iage]=0.0; + + for (i=1; i<=imx; i++) { /* Each individual */ + bool=1; + if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ + for (z1=1; z1<=cptcoveff; z1++) + if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) + bool=0; + } + if (bool==1) { + /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ + for(mi=1; mi=firstpass && m <=lastpass){ y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */ if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ if(agev[m][i]==0) agev[m][i]=iagemax+1; if(agev[m][i]==1) agev[m][i]=iagemax+2; if((int)agev[m][i] iagemax+3) printf("Error on individual =%d agev[m][i]=%f m=%d\n",i, agev[m][i],m); - if (s[m][i]>0 && s[m][i]<=nlstate) { + if (s[m][i]>0 && s[m][i]<=nlstate) { /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/ - prop[s[m][i]][(int)agev[m][i]] += weight[i]; - prop[s[m][i]][iagemax+3] += weight[i]; - } - } + prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */ + prop[s[m][i]][iagemax+3] += weight[i]; + } /* end valid statuses */ + } /* end selection of dates */ } /* end selection of waves */ - } - } - for(i=iagemin; i <= iagemax+3; i++){ - for(jk=1,posprop=0; jk <=nlstate ; jk++) { - posprop += prop[jk][i]; - } - - for(jk=1; jk <=nlstate ; jk++){ - if( i <= iagemax){ - if(posprop>=1.e-5){ - probs[i][jk][j1]= prop[jk][i]/posprop; - } else{ - if(first==1){ - first=0; - printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]); - } + } /* end effective waves */ + } /* end bool */ + } + for(i=iagemin; i <= iagemax+3; i++){ + for(jk=1,posprop=0; jk <=nlstate ; jk++) { + posprop += prop[jk][i]; + } + + for(jk=1; jk <=nlstate ; jk++){ + if( i <= iagemax){ + if(posprop>=1.e-5){ + probs[i][jk][j1]= prop[jk][i]/posprop; + } else{ + if(first==1){ + first=0; + printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]); } - } - }/* end jk */ - }/* end i */ + } + } + }/* end jk */ + }/* end i */ /*} *//* end i1 */ } /* end j1 */ @@ -3208,17 +3598,18 @@ void concatwav(int wav[], int **dh, int int i, mi, m; /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; double sum=0., jmean=0.;*/ - int first; + int first, firstwo; int j, k=0,jk, ju, jl; double sum=0.; first=0; + firstwo=0; jmin=100000; jmax=-1; jmean=0.; - for(i=1; i<=imx; i++){ + for(i=1; i<=imx; i++){ /* For simple cases and if state is death */ mi=0; m=firstpass; - while(s[m][i] <= nlstate){ + while(s[m][i] <= nlstate){ /* a live state */ if(s[m][i]>=1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5) mw[++mi][i]=m; if(m >=lastpass) @@ -3226,13 +3617,26 @@ void concatwav(int wav[], int **dh, int else m++; }/* end while */ - if (s[m][i] > nlstate){ + if (s[m][i] > nlstate){ /* In a death state */ mi++; /* Death is another wave */ /* if(mi==0) never been interviewed correctly before death */ /* Only death is a correct wave */ mw[mi][i]=m; + }else if (andc[i] != 9999) { /* A death occured after lastpass */ + m++; + mi++; + s[m][i]=nlstate+1; /* We are setting the status to the last of non live state */ + mw[mi][i]=m; + nbwarn++; + if(firstwo==0){ + printf("Warning! Death for individual %ld line=%d occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\nOthers in log file only\n",num[i],i,lastpass,nlstate+1, m); + fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\n",num[i],i,lastpass,nlstate+1, m); + firstwo=1; + } + if(firstwo==1){ + fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\n",num[i],i,lastpass,nlstate+1, m); + } } - wav[i]=mi; if(mi==0){ nbwarn++; @@ -3245,7 +3649,9 @@ void concatwav(int wav[], int **dh, int } } /* end mi==0 */ } /* End individuals */ + /* wav and mw are no more changed */ + for(i=1; i<=imx; i++){ for(mi=1; mi File (multiple files are possible if covariates are present): %s\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev)); - fprintf(fichtm,"\n
Probability is computed over estepm=%d months.

\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); + fprintf(fichtm,"\n
Probability is computed over estepm=%d months.

\n", estepm,subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit); /* fprintf(fichtm,"\n
Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year

\n", stepm,YEARM,digitp,digit); */ /* fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.svg\";replot;",digitp,optionfilefiname,digit); */ - fprintf(ficgp,"\nset out;\nset out \"%s%s.svg\";replot;set out;\n",subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); + fprintf(ficgp,"\nset out;\nset out \"%s%s.svg\";replot;set out;\n",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit); free_vector(xp,1,npar); free_matrix(doldm,1,nlstate,1,nlstate); @@ -4097,9 +4507,9 @@ 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, char strstart[]) + 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 *ncvyearp, int ij, char strstart[]) { - /* Variance of prevalence limit */ + /* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ double **dnewm,**doldm; @@ -4107,6 +4517,7 @@ void varprevlim(char fileres[], double * double *xp; double *gp, *gm; double **gradg, **trgradg; + double **mgm, **mgp; double age,agelim; int theta; @@ -4129,6 +4540,8 @@ void varprevlim(char fileres[], double * if (stepm >= YEARM) hstepm=1; nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ gradg=matrix(1,npar,1,nlstate); + mgp=matrix(1,npar,1,nlstate); + mgm=matrix(1,npar,1,nlstate); gp=vector(1,nlstate); gm=vector(1,nlstate); @@ -4136,18 +4549,27 @@ void varprevlim(char fileres[], double * for(i=1; i<=npar; i++){ /* Computes gradient */ xp[i] = x[i] + (i==theta ?delti[theta]:0); } - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); - for(i=1;i<=nlstate;i++) + if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + else + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + for(i=1;i<=nlstate;i++){ gp[i] = prlim[i][i]; - + mgp[theta][i] = prlim[i][i]; + } for(i=1; i<=npar; i++) /* Computes gradient */ xp[i] = x[i] - (i==theta ?delti[theta]:0); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); - for(i=1;i<=nlstate;i++) + if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + else + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + for(i=1;i<=nlstate;i++){ gm[i] = prlim[i][i]; - + mgm[theta][i] = prlim[i][i]; + } for(i=1;i<=nlstate;i++) gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta]; + /* gradg[theta][2]= -gradg[theta][1]; */ /* For testing if nlstate=2 */ } /* End theta */ trgradg =matrix(1,nlstate,1,npar); @@ -4155,11 +4577,34 @@ void varprevlim(char fileres[], double * for(j=1; j<=nlstate;j++) for(theta=1; theta <=npar; theta++) trgradg[j][theta]=gradg[theta][j]; + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */ + /* printf("\nmgm mgp %d ",(int)age); */ + /* for(j=1; j<=nlstate;j++){ */ + /* printf(" %d ",j); */ + /* for(theta=1; theta <=npar; theta++) */ + /* printf(" %d %lf %lf",theta,mgm[theta][j],mgp[theta][j]); */ + /* printf("\n "); */ + /* } */ + /* } */ + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */ + /* printf("\n gradg %d ",(int)age); */ + /* for(j=1; j<=nlstate;j++){ */ + /* printf("%d ",j); */ + /* for(theta=1; theta <=npar; theta++) */ + /* printf("%d %lf ",theta,gradg[theta][j]); */ + /* printf("\n "); */ + /* } */ + /* } */ for(i=1;i<=nlstate;i++) varpl[i][(int)age] =0.; + if((int)age==79 ||(int)age== 80 ||(int)age== 81){ matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); + }else{ + matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); + matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); + } for(i=1;i<=nlstate;i++) varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ @@ -4169,6 +4614,8 @@ void varprevlim(char fileres[], double * fprintf(ficresvpl,"\n"); free_vector(gp,1,nlstate); free_vector(gm,1,nlstate); + free_matrix(mgm,1,npar,1,nlstate); + free_matrix(mgp,1,npar,1,nlstate); free_matrix(gradg,1,npar,1,nlstate); free_matrix(trgradg,1,nlstate,1,npar); } /* End age */ @@ -4200,20 +4647,20 @@ void varprob(char optionfilefiname[], do char fileresprobcor[FILENAMELENGTH]; double ***varpij; - strcpy(fileresprob,"prob"); + strcpy(fileresprob,"PROB_"); strcat(fileresprob,fileres); if((ficresprob=fopen(fileresprob,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileresprob); fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob); } - strcpy(fileresprobcov,"probcov"); - strcat(fileresprobcov,fileres); + strcpy(fileresprobcov,"PROBCOV_"); + strcat(fileresprobcov,fileresu); if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileresprobcov); fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov); } - strcpy(fileresprobcor,"probcor"); - strcat(fileresprobcor,fileres); + strcpy(fileresprobcor,"PROBCOR_"); + strcat(fileresprobcor,fileresu); if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileresprobcor); fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor); @@ -4471,13 +4918,13 @@ To be simple, these graphs help to under fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2); fprintf(ficgp,"\nset ter svg size 640, 480"); fprintf(fichtmcov,"\n
Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\ - :\ -%s%d%1d%1d-%1d%1d.svg, ",k1,l1,k2,l2,\ - subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2,\ - subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); - fprintf(fichtmcov,"\n
",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); + :\ +%s_%d%1d%1d-%1d%1d.svg, ",k1,l1,k2,l2,\ + subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2,\ + subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); + fprintf(fichtmcov,"\n
",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); fprintf(fichtmcov,"\n
Correlation at age %d (%.3f),",(int) age, c12); - fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); + fprintf(ficgp,"\nset out \"%s_%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\ @@ -4494,7 +4941,7 @@ To be simple, these graphs help to under }/* if first */ } /* age mod 5 */ } /* end loop age */ - fprintf(ficgp,"\nset out;\nset out \"%s%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); + fprintf(ficgp,"\nset out;\nset out \"%s_%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); first=1; } /*l12 */ } /* k12 */ @@ -4516,33 +4963,38 @@ To be simple, these graphs help to under /******************* Printing html file ***********/ -void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \ +void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \ int lastpass, int stepm, int weightopt, char model[],\ int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ - int popforecast, int estepm ,\ - double jprev1, double mprev1,double anprev1, \ - double jprev2, double mprev2,double anprev2){ + int popforecast, int prevfcast, int estepm , \ + double jprev1, double mprev1,double anprev1, double dateprev1, \ + double jprev2, double mprev2,double anprev2, double dateprev2){ int jj1, k1, i1, cpt; fprintf(fichtm,""); - fprintf(fichtm,"