--- imach/src/imach.c 2016/08/23 16:51:20 1.234 +++ imach/src/imach.c 2016/08/26 15:51:03 1.239 @@ -1,6 +1,23 @@ -/* $Id: imach.c,v 1.234 2016/08/23 16:51:20 brouard Exp $ +/* $Id: imach.c,v 1.239 2016/08/26 15:51:03 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.239 2016/08/26 15:51:03 brouard + Summary: Improvement in Powell output in order to copy and paste + + Author: + + Revision 1.238 2016/08/26 14:23:35 brouard + Summary: Starting tests of 0.99 + + Revision 1.237 2016/08/26 09:20:19 brouard + Summary: to valgrind + + Revision 1.236 2016/08/25 10:50:18 brouard + *** empty log message *** + + Revision 1.235 2016/08/25 06:59:23 brouard + *** empty log message *** + Revision 1.234 2016/08/23 16:51:20 brouard *** empty log message *** @@ -902,12 +919,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.234 2016/08/23 16:51:20 brouard Exp $ */ +/* $Id: imach.c,v 1.239 2016/08/26 15:51:03 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018"; -char fullversion[]="$Revision: 1.234 $ $Date: 2016/08/23 16:51:20 $"; +char fullversion[]="$Revision: 1.239 $ $Date: 2016/08/26 15:51:03 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1105,6 +1122,16 @@ int *TvarsDind; int *TvarsQ; int *TvarsQind; +#define MAXRESULTLINES 10 +int nresult=0; +int TKresult[MAXRESULTLINES]; +int Tresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */ +int Tinvresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */ +int Tvresult[MAXRESULTLINES][NCOVMAX]; /* For dummy variable , variable # (output) */ +double Tqresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , value (output) */ +double Tqinvresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , value (output) */ +int Tvqresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , variable # (output) */ + /* int *TDvar; /\**< TDvar[1]=4, TDvarF[2]=3, TDvar[3]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 *\/ */ int *TvarF; /**< TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ int *TvarFind; /**< TvarFind[1]=6, TvarFind[2]=7, Tvarind[3]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ @@ -1126,6 +1153,8 @@ double *Tvalsel; /**< Selected modality int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */ int *Fixed; /** Fixed[k] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */ int *Dummy; /** Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */ +int *DummyV; /** Dummy[v] 0=dummy (0 1), 1 quantitative */ +int *FixedV; /** FixedV[v] 0 fixed, 1 varying */ int *Tage; int anyvaryingduminmodel=0; /**< Any varying dummy in Model=1 yes, 0 no, to avoid a loop on waves in freq */ int *Tmodelind; /** Tmodelind[Tvaraff[3]]=9 for V1 position,Tvaraff[1]@9={4, 3, 1, 0, 0, 0, 0, 0, 0}, model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/ @@ -1135,11 +1164,10 @@ int *Ndum; /** Freq of modality (tricode /* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */ int **Tvard; int *Tprod;/**< Gives the k position of the k1 product */ +/* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3 */ int *Tposprod; /**< Gives the k1 product from the k position */ -/* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3 - if V2+V1+V1*V4+age*V3+V3*V2 TProd[k1=2]=5 (V3*V2) - Tposprod[k]=k1 , Tposprod[3]=1, Tposprod[5]=2 -*/ + /* if V2+V1+V1*V4+age*V3+V3*V2 TProd[k1=2]=5 (V3*V2) */ + /* Tposprod[k]=k1 , Tposprod[3]=1, Tposprod[5(V3*V2)]=2 (2nd product without age) */ int cptcovprod, *Tvaraff, *invalidvarcomb; double *lsurv, *lpop, *tpop; @@ -2041,7 +2069,7 @@ void powell(double p[], double **xi, int void linmin(double p[], double xi[], int n, double *fret, double (*func)(double []),int *flat); #endif - int i,ibig,j; + int i,ibig,j,jk,k; double del,t,*pt,*ptt,*xit; double directest; double fp,fptt; @@ -2073,13 +2101,49 @@ void powell(double p[], double **xi, int fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); /* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */ for (i=1;i<=n;i++) { - printf(" %d %.12f",i, p[i]); - fprintf(ficlog," %d %.12lf",i, p[i]); fprintf(ficrespow," %.12lf", p[i]); } + fprintf(ficrespow,"\n");fflush(ficrespow); + printf("\n#model= 1 + age "); + fprintf(ficlog,"\n#model= 1 + age "); + if(nagesqr==1){ + printf(" + age*age ",Tvar[j]); + fprintf(ficlog," + age*age ",Tvar[j]); + } + for(j=1;j <=ncovmodel-2;j++){ + if(Typevar[j]==0) { + printf(" + V%d ",Tvar[j]); + fprintf(ficlog," + V%d ",Tvar[j]); + }else if(Typevar[j]==1) { + printf(" + V%d*age ",Tvar[j]); + fprintf(ficlog," + V%d*age ",Tvar[j]); + }else if(Typevar[j]==2) { + printf(" + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + fprintf(ficlog," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); + } + } printf("\n"); +/* printf("12 47.0114589 0.0154322 33.2424412 0.3279905 2.3731903 */ +/* 13 -21.5392400 0.1118147 1.2680506 1.2973408 -1.0663662 */ fprintf(ficlog,"\n"); - fprintf(ficrespow,"\n");fflush(ficrespow); + for(i=1,jk=1; i <=nlstate; i++){ + for(k=1; k <=(nlstate+ndeath); k++){ + if (k != i) { + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + fprintf(ficres,"%1d%1d ",i,k); + for(j=1; j <=ncovmodel; j++){ + printf("%12.7f ",p[jk]); + fprintf(ficlog,"%12.7f ",p[jk]); + fprintf(ficres,"%12.7f ",p[jk]); + jk++; + } + printf("\n"); + fprintf(ficlog,"\n"); + fprintf(ficres,"\n"); + } + } + } if(*iter <=3){ tml = *localtime(&rcurr_time); strcpy(strcurr,asctime(&tml)); @@ -2333,9 +2397,11 @@ 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 *ncvyear, int ij) + double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij, int nres) { - /* Computes the prevalence limit in each live state at age x and for covariate combiation ij by left multiplying the unit + /* Computes the prevalence limit in each live state at age x and for covariate combination ij + (and selected quantitative values in nres) + by left multiplying the unit 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 */ @@ -2387,27 +2453,36 @@ double **prevalim(double **prlim, int nl for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ /* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */ cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; - printf("prevalim ij=%d k=%d TvarsD[%d]=%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); + /* printf("prevalim Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ } for (k=1; k<=nsq;k++) { /* For single varying covariates only */ /* Here comes the value of quantitative after renumbering k with single quantitative covariates */ - /* cov[2+nagesqr+TvarsQind[k]]=qselvar[k]; */ - printf("prevalim ij=%d k=%d TvarsQind[%d]=%d \n",ij,k,k,TvarsQind[k]); + cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; + /* printf("prevalim Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][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]][codtabm(ij,Tvar[k])]*cov[2]; */ - for (k=1; k<=cptcovage;k++){ + for (k=1; k<=cptcovage;k++){ /* For product with age */ if(Dummy[Tvar[Tage[k]]]){ cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; } else{ - ; - /* cov[2+nagesqr+Tage[k]]=qselvar[k]; */ + cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; } - printf("prevalim Age ij=%d k=%d Tage[%d]=%d \n",ij,k,k,Tage[k]); + /* printf("prevalim Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */ } - for (k=1; k<=cptcovprod;k++){ /* */ - printf("prevalim Prod ij=%d k=%d Tprod[%d]=%d Tvard[%d][1]=%d, Tvard[%d][2]=%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; + for (k=1; k<=cptcovprod;k++){ /* For product without age */ + /* printf("prevalim Prod ij=%d k=%d Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */ + if(Dummy[Tvard[k][1]==0]){ + if(Dummy[Tvard[k][2]==0]){ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; + }else{ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; + } + }else{ + if(Dummy[Tvard[k][2]==0]){ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; + }else{ + cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; + } + } } /*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]);*/ @@ -2530,8 +2605,6 @@ Earliest age to start was %d-%d=%d, ncvl 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]][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]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ @@ -2856,7 +2929,7 @@ double **matprod2(double **out, double * /************* Higher Matrix Product ***************/ -double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij ) +double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij, int nres ) { /* Computes the transition matrix starting at age 'age' and combination of covariate values corresponding to ij over 'nhstepm*hstepm*stepm' months (i.e. until @@ -2892,16 +2965,34 @@ double ***hpxij(double ***po, int nhstep cov[2]=agexact; if(nagesqr==1) cov[3]= agexact*agexact; - for (k=1; k<=cptcovn;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]]][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]][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])]; */ + for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ + /* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */ + cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; + /* printf("hpxij Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ + } + for (k=1; k<=nsq;k++) { /* For single varying covariates only */ + /* Here comes the value of quantitative after renumbering k with single quantitative covariates */ + cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; + /* printf("hPxij Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */ + } + for (k=1; k<=cptcovage;k++){ + if(Dummy[Tvar[Tage[k]]]){ + cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; + } else{ + cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; + } + /* printf("hPxij Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */ + } + for (k=1; k<=cptcovprod;k++){ /* */ + /* printf("hPxij Prod ij=%d k=%d Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; + } + /* for (k=1; k<=cptcovn;k++) */ + /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */ + /* for (k=1; k<=cptcovage;k++) /\* Should start at cptcovn+1 *\/ */ + /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */ + /* for (k=1; k<=cptcovprod;k++) /\* Useless because included in cptcovn *\/ */ + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; */ /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ @@ -4058,7 +4149,7 @@ void freqsummary(char fileres[], int ia 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); + fprintf(ficresphtm,"Current page is file %s
\n\n

Frequencies and prevalence by age at begin of transition and dummy covariate value at beginning of transition

\n",fileresphtm, fileresphtm); strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { @@ -4862,7 +4953,7 @@ void concatwav(int wav[], int **dh, int /*********** Health Expectancies ****************/ -void evsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] ) + void evsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[], int nres ) { /* Health expectancies, no variances */ @@ -4872,7 +4963,7 @@ void evsij(double ***eij, double x[], in double ***p3mat; double eip; - pstamp(ficreseij); + /* pstamp(ficreseij); */ fprintf(ficreseij,"# (a) Life expectancies by health status at initial age and (b) health expectancies by health status at initial age\n"); fprintf(ficreseij,"# Age"); for(i=1; i<=nlstate;i++){ @@ -4935,7 +5026,7 @@ void evsij(double ***eij, double x[], in /* Computed by stepm unit matrices, product of hstepma matrices, stored in an array of nhstepma length: nhstepma=10, hstepm=4, stepm=6 months */ - hpxij(p3mat,nhstepma,age,hstepm,x,nlstate,stepm,oldm, savm, cij); + hpxij(p3mat,nhstepma,age,hstepm,x,nlstate,stepm,oldm, savm, cij, nres); hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ @@ -4970,7 +5061,7 @@ void evsij(double ***eij, double x[], in } -void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] ) + void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[], int nres ) { /* Covariances of health expectancies eij and of total life expectancies according @@ -5083,8 +5174,8 @@ void cvevsij(double ***eij, double x[], xp[i] = x[i] + (i==theta ?delti[theta]:0); xm[i] = x[i] - (i==theta ?delti[theta]:0); } - hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij); - hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij); + hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij, nres); + hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij, nres); for(j=1; j<= nlstate; j++){ for(i=1; i<=nlstate; i++){ @@ -5125,7 +5216,7 @@ void cvevsij(double ***eij, double x[], } /* Computing expectancies */ - hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); + hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij,nres); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ @@ -5180,7 +5271,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 *ncvyearp, 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[], int nres) { /* Variance of health expectancies */ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ @@ -5237,6 +5328,14 @@ void cvevsij(double ***eij, double x[], 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); + fprintf(ficresprobmorprev,"# Selected quantitative variables and dummies"); + for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ + fprintf(ficresprobmorprev," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); + } + for(j=1;j<=cptcoveff;j++) + fprintf(ficresprobmorprev,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(ij,j)]); + fprintf(ficresprobmorprev,"\n"); + fprintf(ficresprobmorprev,"# Age cov=%-d",ij); for(j=nlstate+1; j<=(nlstate+ndeath);j++){ fprintf(ficresprobmorprev," p.%-d SE",j); @@ -5306,7 +5405,7 @@ void cvevsij(double ***eij, double x[], xp[i] = x[i] + (i==theta ?delti[theta]:0); } - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nresult); if (popbased==1) { if(mobilav ==0){ @@ -5318,7 +5417,7 @@ void cvevsij(double ***eij, double x[], } } - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); /* Returns p3mat[i][j][h] for h=1 to nhstepm */ + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* 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++) @@ -5338,7 +5437,7 @@ void cvevsij(double ***eij, double x[], for(i=1; i<=npar; i++) /* Computes gradient x - delta */ xp[i] = x[i] - (i==theta ?delti[theta]:0); - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij); + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nresult); if (popbased==1) { if(mobilav ==0){ @@ -5350,7 +5449,7 @@ void cvevsij(double ***eij, double x[], } } - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */ for(h=0; h<=nhstepm; h++){ @@ -5415,7 +5514,7 @@ void cvevsij(double ***eij, double x[], /* end ppptj */ /* x centered again */ - prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij); + prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij, nresult); if (popbased==1) { if(mobilav ==0){ @@ -5431,7 +5530,7 @@ void cvevsij(double ***eij, double x[], 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); + hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij, nres); 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]; @@ -5494,7 +5593,7 @@ void cvevsij(double ***eij, double x[], } /* 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 *ncvyearp, 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[], int nres) { /* 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);*/ @@ -5537,9 +5636,9 @@ void cvevsij(double ***eij, double x[], xp[i] = x[i] + (i==theta ?delti[theta]:0); } if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); else - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); for(i=1;i<=nlstate;i++){ gp[i] = prlim[i][i]; mgp[theta][i] = prlim[i][i]; @@ -5547,9 +5646,9 @@ void cvevsij(double ***eij, double x[], for(i=1; i<=npar; i++) /* Computes gradient */ xp[i] = x[i] - (i==theta ?delti[theta]:0); if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); else - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); for(i=1;i<=nlstate;i++){ gm[i] = prlim[i][i]; mgm[theta][i] = prlim[i][i]; @@ -5957,11 +6056,13 @@ void printinghtml(char fileresu[], char int popforecast, int prevfcast, int backcast, int estepm , \ double jprev1, double mprev1,double anprev1, double dateprev1, \ double jprev2, double mprev2,double anprev2, double dateprev2){ - int jj1, k1, i1, cpt; + int jj1, k1, i1, cpt, k4, nres; fprintf(fichtm,""); + fprintf(fichtm,"", model); fprintf(fichtm,"