--- imach/src/imach.c 2015/08/18 13:32:00 1.194 +++ imach/src/imach.c 2015/11/17 22:12:03 1.209 @@ -1,6 +1,72 @@ -/* $Id: imach.c,v 1.194 2015/08/18 13:32:00 brouard Exp $ +/* $Id: imach.c,v 1.209 2015/11/17 22:12:03 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + 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 + + Even model=1+age+V2. did not work anymore + + Revision 1.199 2015/09/07 14:09:23 brouard + Summary: 0.98q6 changing default small png format for graph to vectorized svg. + + Revision 1.198 2015/09/03 07:14:39 brouard + Summary: 0.98q5 Flavia + + Revision 1.197 2015/09/01 18:24:39 brouard + *** empty log message *** + + Revision 1.196 2015/08/18 23:17:52 brouard + Summary: 0.98q5 + + Revision 1.195 2015/08/18 16:28:39 brouard + Summary: Adding a hack for testing purpose + + After reading the title, ftol and model lines, if the comment line has + a q, starting with #q, the answer at the end of the run is quit. It + permits to run test files in batch with ctest. The former workaround was + $ echo q | imach foo.imach + Revision 1.194 2015/08/18 13:32:00 brouard Summary: Adding error when the covariance matrix doesn't contain the exact number of lines required by the model line. @@ -613,6 +679,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 *\/ */ @@ -683,7 +753,7 @@ typedef struct { #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */ #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) ; +#define codtabm(h,k) (1 & (h-1) >> (k-1))+1 #define MAXN 20000 #define YEARM 12. /**< Number of months per year */ #define AGESUP 130 @@ -700,11 +770,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.194 2015/08/18 13:32:00 brouard Exp $ */ +/* $Id: imach.c,v 1.209 2015/11/17 22:12:03 brouard Exp $ */ /* $State: Exp $ */ - -char version[]="Imach version 0.98q5, August 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.194 $ $Date: 2015/08/18 13:32:00 $"; +#include "version.h" +char version[]=__IMACH_VERSION__; +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.209 $ $Date: 2015/11/17 22:12:03 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -770,7 +841,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]; @@ -838,7 +909,7 @@ int estepm; int m,nb; long *num; -int firstpass=0, lastpass=4,*cod, *Tage,*cens; +int firstpass=0, lastpass=4,*cod, *cens; int *ncodemax; /* ncodemax[j]= Number of modalities of the j th covariate for which somebody answered excluding undefined. Usually 2: 0 and 1. */ @@ -858,8 +929,9 @@ double **covar; /**< covar[j,i], value * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */ double idx; int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ +int *Tage; int *Ndum; /** Freq of modality (tricode */ -int **codtab; /**< codtab=imatrix(1,100,1,10); */ +/* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */ int **Tvard, *Tprod, cptcovprod, *Tvaraff; double *lsurv, *lpop, *tpop; @@ -893,7 +965,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 ); @@ -1565,23 +1637,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)) */ @@ -1589,14 +1667,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]) */ @@ -1609,22 +1695,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); @@ -1658,7 +1759,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; @@ -1702,10 +1803,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 */ @@ -1803,7 +1904,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); @@ -1818,9 +1919,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 */ @@ -1828,17 +1929,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++) { @@ -1868,18 +1975,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); @@ -1888,20 +2017,25 @@ 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; if(nagesqr==1) cov[3]= agefin*agefin;; for (k=1; k<=cptcovn;k++) { - cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; - /*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/ + /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; + /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */ } /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ - for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]*cov[2]; + /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */ + for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2]; for (k=1; k<=cptcovprod;k++) /* Useless */ - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/ @@ -1912,25 +2046,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 */ } @@ -2074,12 +2228,15 @@ double ***hpxij(double ***po, int nhstep if(nagesqr==1) cov[3]= agexact*agexact; for (k=1; k<=cptcovn;k++) - cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; + /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */ /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ - cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; + cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; + /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */ for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */ - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ @@ -2503,9 +2660,9 @@ 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,"%9ld %6d %2d %2d %1d %1d %3d %11.6f %8.4f\ + fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\ %11.6f %11.6f %11.6f ", \ - num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i], + num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); for(k=1,llt=0.,l=0.; k<=nlstate; k++){ llt +=ll[k]*gipmx/gsw; @@ -2537,14 +2694,14 @@ void likelione(FILE *ficres,double p[], int k; if(*globpri !=0){ /* Just counts and sums, no printings */ - strcpy(fileresilk,"ilk"); - strcat(fileresilk,fileres); + strcpy(fileresilk,"ILK_"); + strcat(fileresilk,fileresu); if((ficresilk=fopen(fileresilk,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileresilk); fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk); } - fprintf(ficresilk, "#individual(line's_record) s1 s2 wave# effective_wave# number_of_matrices_product pij weight -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n"); - fprintf(ficresilk, "#num_i i s1 s2 mi mw dh likeli weight 2wlli out sav "); + fprintf(ficresilk, "#individual(line's_record) count age s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n"); + fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav "); /* 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; k<=nlstate; k++) fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k); @@ -2554,9 +2711,23 @@ void likelione(FILE *ficres,double p[], *fretone=(*funcone)(p); if(*globpri !=0){ fclose(ficresilk); - fprintf(fichtm,"\n
File 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 %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; } @@ -2587,7 +2758,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); @@ -2630,32 +2801,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); @@ -2666,9 +2837,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]);*/ @@ -2702,53 +2873,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; @@ -2764,9 +2960,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 @@ -2780,48 +2976,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) { @@ -2898,7 +3171,7 @@ 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); + fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart); } /************ Frequencies ********************/ @@ -2914,8 +3187,8 @@ void freqsummary(char fileres[], int ia pp=vector(1,nlstate); prop=matrix(1,nlstate,iagemin,iagemax+3); - strcpy(fileresp,"p"); - strcat(fileresp,fileres); + strcpy(fileresp,"P_"); + strcat(fileresp,fileresu); if((ficresp=fopen(fileresp,"w"))==NULL) { printf("Problem with prevalence resultfile: %s\n", fileresp); fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); @@ -2950,13 +3223,13 @@ void freqsummary(char fileres[], int ia 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]][codtab[j1][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 */ bool=0; - /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtab[%d][%d]=%d, nbcode[Tvaraff][codtab[%d][%d]=%d, j1=%d\n", - bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtab[j1][z1], - j1,z1,nbcode[Tvaraff[z1]][codtab[j1][z1]],j1);*/ - /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtab[7][3]=1 and nbcde[3][?]=1*/ + /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", + bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), + 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*/ } } @@ -2985,10 +3258,10 @@ void freqsummary(char fileres[], int ia 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]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficresp, "**********\n#"); fprintf(ficlog, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficlog, "**********\n#"); } for(i=1; i<=nlstate;i++) @@ -3116,7 +3389,7 @@ void prevalence(double ***probs, double bool=1; if (cptcovn>0) { for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) + if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) bool=0; } if (bool==1) { @@ -3382,14 +3655,15 @@ void tricode(int *Tvar, int **nbcode, in nbcode[Tvar[j]][1]=0; nbcode[Tvar[j]][2]=1; nbcode[Tvar[j]][3]=2; + To be continued (not working yet). */ - ij=0; /* ij is similar to i but can jumps over null modalities */ - for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 to 1*/ - if (Ndum[i] == 0) { /* If at least one individual responded to this modality k */ + ij=0; /* ij is similar to i but can jump over null modalities */ + for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ + if (Ndum[i] == 0) { /* If nobody responded to this modality k */ break; } ij++; - nbcode[Tvar[j]][ij]=i; /* stores the original modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/ + nbcode[Tvar[j]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/ cptcode = ij; /* New max modality for covar j */ } /* end of loop on modality i=-1 to 1 or more */ @@ -3753,7 +4027,7 @@ void cvevsij(double ***eij, double x[], } /************ Variance ******************/ -void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[]) + void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[]) { /* Variance of health expectancies */ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ @@ -3782,11 +4056,11 @@ void varevsij(char optionfilefiname[], d if(popbased==1){ if(mobilav!=0) - strcpy(digitp,"-populbased-mobilav-"); - else strcpy(digitp,"-populbased-nomobil-"); + strcpy(digitp,"-POPULBASED-MOBILAV_"); + else strcpy(digitp,"-POPULBASED-NOMOBIL_"); } else - strcpy(digitp,"-stablbased-"); + strcpy(digitp,"-STABLBASED_"); if (mobilav!=0) { mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); @@ -3796,18 +4070,17 @@ void varevsij(char optionfilefiname[], d } } - strcpy(fileresprobmorprev,"prmorprev"); + strcpy(fileresprobmorprev,"PRMORPREV-"); sprintf(digit,"%-d",ij); /*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/ strcat(fileresprobmorprev,digit); /* Tvar to be done */ strcat(fileresprobmorprev,digitp); /* Popbased or not, mobilav or not */ - strcat(fileresprobmorprev,fileres); + strcat(fileresprobmorprev,fileresu); if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileresprobmorprev); fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev); } printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); - fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); pstamp(ficresprobmorprev); fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm); @@ -3818,8 +4091,10 @@ void varevsij(char optionfilefiname[], d fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j); } fprintf(ficresprobmorprev,"\n"); + fprintf(ficgp,"\n# Routine varevsij"); - /* fprintf(fichtm, "#Local time at start: %s", strstart);*/ + fprintf(ficgp,"\nunset title \n"); +/* fprintf(fichtm, "#Local time at start: %s", strstart);*/ fprintf(fichtm,"\n

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

    • \n"); fprintf(fichtm,"\n
    %s
    \n",digitp); /* } */ @@ -3854,9 +4129,9 @@ void varevsij(char optionfilefiname[], d /* 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 function hpijx to understand why (it is linked to memory size questions) */ - /* We decided (b) to get a life expectancy respecting the most precise curvature of the + nstepm is the number of stepm from age to agelim. + Look at function hpijx to understand why because of memory size limitations, + 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 every two years of age and if you sum them up and add 1 year (area under the trapezoids) you won't get the same @@ -3877,8 +4152,8 @@ void varevsij(char optionfilefiname[], d for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/ xp[i] = x[i] + (i==theta ?delti[theta]:0); } - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); + + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); if (popbased==1) { if(mobilav ==0){ @@ -3890,13 +4165,14 @@ void varevsij(char optionfilefiname[], d } } + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); /* Returns p3mat[i][j][h] for h=1 to nhstepm */ for(j=1; j<= nlstate; j++){ for(h=0; h<=nhstepm; h++){ for(i=1, gp[h][j]=0.;i<=nlstate;i++) gp[h][j] += prlim[i][i]*p3mat[i][j][h]; } } - /* This for computing probability of death (h=1 means + /* Next for computing probability of death (h=1 means computed over hstepm matrices product = hstepm*stepm months) as a weighted average of prlim. */ @@ -3908,8 +4184,8 @@ void varevsij(char optionfilefiname[], d for(i=1; i<=npar; i++) /* Computes gradient x - delta */ xp[i] = x[i] - (i==theta ?delti[theta]:0); - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); + + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij); if (popbased==1) { if(mobilav ==0){ @@ -3921,6 +4197,8 @@ void varevsij(char optionfilefiname[], d } } + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); + for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */ for(h=0; h<=nhstepm; h++){ for(i=1, gm[h][j]=0.;i<=nlstate;i++) @@ -3983,8 +4261,8 @@ void varevsij(char optionfilefiname[], d varppt[j][i]=doldmp[j][i]; /* end ppptj */ /* x centered again */ - hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); - prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij); + + prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij); if (popbased==1) { if(mobilav ==0){ @@ -4000,6 +4278,7 @@ void varevsij(char optionfilefiname[], d computed over hstepm (estepm) matrices product = hstepm*stepm months) as a weighted average of prlim. */ + hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); for(j=nlstate+1;j<=nlstate+ndeath;j++){ for(i=1,gmp[j]=0.;i<= nlstate; i++) gmp[j] += prlim[i][i]*p3mat[i][j][1]; @@ -4031,9 +4310,11 @@ void varevsij(char optionfilefiname[], d free_vector(gmp,nlstate+1,nlstate+ndeath); free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath); free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/ - fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); + /* fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); */ + fprintf(ficgp,"\nunset parametric;unset label; set ter svg size 640, 480"); /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */ fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); + fprintf(ficgp,"\nset out \"%s%s.svg\";",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit); /* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ @@ -4041,11 +4322,11 @@ void varevsij(char optionfilefiname[], d fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev)); fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev)); fprintf(fichtm,"\n
    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 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(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.png\";replot;",digitp,optionfilefiname,digit); */ - fprintf(ficgp,"\nset out \"%s%s.png\";replot;\n",subdirf3(optionfilefiname,"varmuptjgr",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); free_vector(xp,1,npar); free_matrix(doldm,1,nlstate,1,nlstate); @@ -4060,9 +4341,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; @@ -4070,6 +4351,7 @@ void varprevlim(char fileres[], double * double *xp; double *gp, *gm; double **gradg, **trgradg; + double **mgm, **mgp; double age,agelim; int theta; @@ -4092,6 +4374,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); @@ -4099,18 +4383,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); @@ -4118,11 +4411,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 */ @@ -4132,6 +4448,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 */ @@ -4163,20 +4481,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); @@ -4218,10 +4536,9 @@ void varprob(char optionfilefiname[], do fprintf(fichtm,"\n

  • Computing and drawing one step probabilities with their confidence intervals

    • \n"); fprintf(fichtm,"\n"); - fprintf(fichtm,"\n

  • Matrix of variance-covariance of pairs of step probabilities (drawings)

    • \n",optionfilehtmcov); - fprintf(fichtmcov,"\n

    Matrix of variance-covariance of pairs of step probabilities

    \n\ - file %s
    \n",optionfilehtmcov); - fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (pij, pkl) are estimated\ + fprintf(fichtm,"\n

  • Matrix of variance-covariance of one-step probabilities (drawings)

    this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.
    • \n",optionfilehtmcov); + fprintf(fichtmcov,"Current page is file %s
    \n\n

    Matrix of variance-covariance of pairs of step probabilities

    \n",optionfilehtmcov, optionfilehtmcov); + fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (pij, pkl) are estimated \ and drawn. It helps understanding how is the covariance between two incidences.\ They are expressed in year-1 in order to be less dependent of stepm.
    \n"); fprintf(fichtmcov,"\n
    Contour plot corresponding to x'cov-1x = 4 (where x is the column vector (pij,pkl)) are drawn. \ @@ -4242,23 +4559,23 @@ To be simple, these graphs help to under /*j1++;*/ if (cptcovn>0) { fprintf(ficresprob, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficresprob, "**********\n#\n"); fprintf(ficresprobcov, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficresprobcov, "**********\n#\n"); fprintf(ficgp, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficgp, "**********\n#\n"); fprintf(fichtmcov, "\n
    ********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(fichtmcov, "**********\n
    "); fprintf(ficresprobcor, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); fprintf(ficresprobcor, "**********\n#"); } @@ -4271,7 +4588,8 @@ To be simple, these graphs help to under if(nagesqr==1) cov[3]= age*age; for (k=1; k<=cptcovn;k++) { - cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];/* j1 1 2 3 4 + cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; + /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4 * 1 1 1 1 1 * 2 2 1 1 1 * 3 1 2 1 1 @@ -4279,9 +4597,9 @@ To be simple, these graphs help to under /* nbcode[1][1]=0 nbcode[1][2]=1;*/ } /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ - for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; + for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; for (k=1; k<=cptcovprod;k++) - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; for(theta=1; theta <=npar; theta++){ @@ -4429,17 +4747,18 @@ To be simple, these graphs help to under /* mu2+ v21*lc1*cost + v22*lc2*sin(t) */ if(first==1){ first=0; + fprintf(ficgp,"\n# Ellipsoids of confidence\n#\n"); fprintf(ficgp,"\nset parametric;unset label"); 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 png small size 320, 240"); + 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.png, ",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.png\"",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",\ @@ -4456,7 +4775,7 @@ To be simple, these graphs help to under }/* if first */ } /* age mod 5 */ } /* end loop age */ - fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\";replot;",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 */ @@ -4478,7 +4797,7 @@ 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 ,\ @@ -4491,20 +4810,20 @@ void printinghtml(char fileres[], char t "); fprintf(fichtm,"