--- imach/src/imach.c 2022/07/25 14:27:23 1.325 +++ imach/src/imach.c 2022/08/07 05:40:09 1.331 @@ -1,6 +1,26 @@ -/* $Id: imach.c,v 1.325 2022/07/25 14:27:23 brouard Exp $ +/* $Id: imach.c,v 1.331 2022/08/07 05:40:09 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.331 2022/08/07 05:40:09 brouard + *** empty log message *** + + Revision 1.330 2022/08/06 07:18:25 brouard + Summary: last 0.99r31 + + * imach.c (Module): Version of imach using partly decoderesult to rebuild xpxij function + + Revision 1.329 2022/08/03 17:29:54 brouard + * imach.c (Module): Many errors in graphs fixed with Vn*age covariates. + + Revision 1.328 2022/07/27 17:40:48 brouard + Summary: valgrind bug fixed by initializing to zero DummyV as well as Tage + + Revision 1.327 2022/07/27 14:47:35 brouard + Summary: Still a problem for one-step probabilities in case of quantitative variables + + Revision 1.326 2022/07/26 17:33:55 brouard + Summary: some test with nres=1 + Revision 1.325 2022/07/25 14:27:23 brouard Summary: r30 @@ -1183,7 +1203,7 @@ typedef struct { #define GNUPLOTPROGRAM "gnuplot" /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ -#define FILENAMELENGTH 132 +#define FILENAMELENGTH 256 #define GLOCK_ERROR_NOPATH -1 /* empty path */ #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ @@ -1218,24 +1238,25 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.325 2022/07/25 14:27:23 brouard Exp $ */ +/* $Id: imach.c,v 1.331 2022/08/07 05:40:09 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; char copyright[]="July 2022,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022"; -char fullversion[]="$Revision: 1.325 $ $Date: 2022/07/25 14:27:23 $"; +char fullversion[]="$Revision: 1.331 $ $Date: 2022/08/07 05:40:09 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ int nagesqr=0, nforce=0; /* nagesqr=1 if model is including age*age, number of forces */ -/* Number of covariates model=V2+V1+ V3*age+V2*V4 */ -int cptcovn=0; /**< cptcovn number of covariates added in the model (excepting constant and age and age*product) */ -int cptcovt=0; /**< cptcovt number of covariates added in the model (excepting constant and age) */ -int cptcovs=0; /**< cptcovs number of simple covariates in the model V2+V1 =2 */ +/* Number of covariates model (1)=V2+V1+ V3*age+V2*V4 */ +/* Model(2) V1 + V2 + V3 + V8 + V7*V8 + V5*V6 + V8*age + V3*age + age*age */ +int cptcovn=0; /**< cptcovn decodemodel: number of covariates k of the models excluding age*products =6 and age*age */ +int cptcovt=0; /**< cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */ +int cptcovs=0; /**< cptcovs number of simple covariates in the model V2+V1 =2 (dummy or quantit or time varying) */ int cptcovsnq=0; /**< cptcovsnq number of simple covariates in the model but non quantitative V2+V1 =2 */ int cptcovage=0; /**< Number of covariates with age: V3*age only =1 */ int cptcovprodnoage=0; /**< Number of covariate products without age */ -int cptcoveff=0; /* Total number of covariates to vary for printing results */ +int cptcoveff=0; /* Total number of covariates to vary for printing results (2**cptcoveff combinations of dummies)(computed in tricode as cptcov) */ int ncovf=0; /* Total number of effective fixed covariates (dummy or quantitative) in the model */ int ncovv=0; /* Total number of effective (wave) varying covariates (dummy or quantitative) in the model */ int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (dummy of quantitative) in the model */ @@ -1409,7 +1430,7 @@ int **nbcode, *Tvar; /**< model=V2 => Tv * V1 V2 V3 V4 V5 V6 V7 V8 Weight ddb ddth d1st s1 V9 V10 V11 V12 s2 V9 V10 V11 V12 * < ncovcol=6 > nqv=2 (V7 V8) dv dv dv qtv dv dv dvv qtv * ntv=3 nqtv=1 - * cptcovn number of covariates (not including constant and age) = # of + plus 1 = 10+1=11 + * cptcovn number of covariates (not including constant and age or age*age) = number of plus sign + 1 = 10+1=11 * For time varying covariate, quanti or dummies * cotqvar[wav][iv(1 to nqtv)][i]= [1][12][i]=(V12) quanti * cotvar[wav][ntv+iv][i]= [3+(1 to nqtv)][i]=(V12) quanti @@ -1435,6 +1456,7 @@ int **nbcode, *Tvar; /**< model=V2 => Tv /* with age product, 3 quant with age product*/ /*Tvar[k]= 5 4 3 6 5 2 7 1 1 */ /* nsd 1 2 3 */ /* Counting single dummies covar fixed or tv */ +/*TnsdVar[Tvar] 1 2 3 */ /*TvarsD[nsd] 4 3 1 */ /* ID of single dummy cova fixed or timevary*/ /*TvarsDind[k] 2 3 9 */ /* position K of single dummy cova */ /* nsq 1 2 */ /* Counting single quantit tv */ @@ -1444,6 +1466,7 @@ int **nbcode, *Tvar; /**< model=V2 => Tv /* cptcovage 1 2 */ /* Counting cov*age in the model equation */ /* Tage[cptcovage]=k 5 8 */ /* Position in the model of ith cov*age */ /* Tvard[1][1]@4={4,3,1,2} V4*V3 V1*V2 */ /* Position in model of the ith prod without age */ +/* Tvardk[4][1]=4;Tvardk[4][2]=3;Tvardk[7][1]=1;Tvardk[7][2]=2 */ /* Variables of a prod at position in the model equation*/ /* TvarF TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 ID of fixed covariates or product V2, V1*V2, V1 */ /* TvarFind; TvarFind[1]=6, TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod) */ /* Type */ @@ -1452,6 +1475,7 @@ int **nbcode, *Tvar; /**< model=V2 => Tv /* D Q D D Q */ /* */ int *TvarsD; +int *TnsdVar; int *TvarsDind; int *TvarsQ; int *TvarsQind; @@ -1460,8 +1484,10 @@ int *TvarsQind; int nresult=0; int parameterline=0; /* # of the parameter (type) line */ int TKresult[MAXRESULTLINESPONE]; +int resultmodel[MAXRESULTLINESPONE][NCOVMAX];/* resultmodel[k1]=k3: k1th position in the model correspond to the k3 position in the resultline */ int Tresult[MAXRESULTLINESPONE][NCOVMAX];/* For dummy variable , value (output) */ int Tinvresult[MAXRESULTLINESPONE][NCOVMAX];/* For dummy variable , value (output) */ +int TinvDoQresult[MAXRESULTLINESPONE][NCOVMAX];/* For dummy variable or quanti value (output) */ int Tvresult[MAXRESULTLINESPONE][NCOVMAX]; /* For dummy variable , variable # (output) */ double Tqresult[MAXRESULTLINESPONE][NCOVMAX]; /* For quantitative variable , value (output) */ double Tqinvresult[MAXRESULTLINESPONE][NCOVMAX]; /* For quantitative variable , value (output) */ @@ -1502,6 +1528,7 @@ int *TmodelInvQind; /** Tmodelqind[1]=1 int *Ndum; /** Freq of modality (tricode */ /* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */ int **Tvard; +int **Tvardk; 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 */ @@ -2793,7 +2820,7 @@ void powell(double p[], double **xi, int } 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)]; + cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,TvarsD[k])]; /* cov[++k1]=nbcode[TvarsD[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)); */ } @@ -2805,7 +2832,7 @@ void powell(double p[], double **xi, int } for (k=1; k<=cptcovage;k++){ /* For product with age */ if(Dummy[Tage[k]]==2){ /* dummy with age */ - 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]; /* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */ } else if(Dummy[Tage[k]]==3){ /* quantitative with age */ cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; @@ -2815,17 +2842,17 @@ void powell(double p[], double **xi, int } 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)]; + if(Dummy[Tvard[k][1]]==0){ + if(Dummy[Tvard[k][2]]==0){ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; /* cov[++k1]=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]; + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * Tqresult[nres][k]; /* cov[++k1]=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]]; + if(Dummy[Tvard[k][2]]==0){ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])] * Tqinvresult[nres][Tvard[k][1]]; /* cov[++k1]=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]]; @@ -2966,7 +2993,7 @@ void powell(double p[], double **xi, int } 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)]; + cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,TvarsD[k])]; /* printf("bprevalim Dummy agefin=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agefin,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ } /* for (k=1; k<=cptcovn;k++) { */ @@ -2986,7 +3013,7 @@ void powell(double p[], double **xi, int for (k=1; k<=cptcovage;k++){ /* For product with age */ /* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age *\/ ERROR ???*/ if(Dummy[Tage[k]]== 2){ /* dummy with age */ - 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]; } else if(Dummy[Tage[k]]== 3){ /* quantitative with age */ cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; } @@ -2994,15 +3021,15 @@ void powell(double p[], double **xi, int } 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)]; + if(Dummy[Tvard[k][1]]==0){ + if(Dummy[Tvard[k][2]]==0){ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; }else{ - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * 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]]; + if(Dummy[Tvard[k][2]]==0){ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])] * Tqinvresult[nres][Tvard[k][1]]; }else{ cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; } @@ -3118,7 +3145,7 @@ double **pmij(double **ps, double *cov, /* printf("Int ji s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */ } ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */ + /* printf("Debug pmij() i=%d j=%d nc=%d s1=%.17f, lnpijopii=%.17f\n",i,j,nc, s1,lnpijopii); */ } } @@ -3134,11 +3162,11 @@ double **pmij(double **ps, double *cov, s1=0; for(j=1; ji} pij/pii=(1-pii)/pii and thus pii is known from s1 */ ps[i][i]=1./(s1+1.); @@ -3394,7 +3422,7 @@ double ***hpxij(double ***po, int nhstep */ - int i, j, d, h, k; + int i, j, d, h, k, k1; double **out, cov[NCOVMAX+1]; double **newm; double agexact; @@ -3417,8 +3445,13 @@ double ***hpxij(double ***po, int nhstep if(nagesqr==1){ cov[3]= agexact*agexact; } - 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 */ + /* Model(2) V1 + V2 + V3 + V8 + V7*V8 + V5*V6 + V8*age + V3*age + age*age */ + /* total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age */ + for(k1=1;k1<=cptcovt;k1++){ /* loop on model equation (including products) */ + if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Single dummy */ +/* V(Tvarsel)=Tvalsel=Tresult[nres][pos](value); V(Tvresult[nres][pos] (variable): V(variable)=value) */ +/* 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 *\/ */ /* codtabm(ij,k) (1 & (ij-1) >> (k-1))+1 */ /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* k 1 2 3 4 5 6 7 8 9 */ @@ -3426,42 +3459,78 @@ double ***hpxij(double ***po, int nhstep /* nsd 1 2 3 */ /* Counting single dummies covar fixed or tv */ /*TvarsD[nsd] 4 3 1 */ /* ID of single dummy cova fixed or timevary*/ /*TvarsDind[k] 2 3 9 */ /* position K of single dummy cova */ - 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++){ /* For product with age V1+V1*age +V4 +age*V3 */ - /* 1+2 Tage[1]=2 TVar[2]=1 Dummy[2]=2, Tage[2]=4 TVar[4]=3 Dummy[4]=3 quant*/ - /* */ - if(Dummy[Tage[k]]== 2){ /* dummy with age */ - /* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age *\/ */ - cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; - } else if(Dummy[Tage[k]]== 3){ /* quantitative with age */ - 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++){ /* For product without age */ - /* 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)]; */ - 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]]; - } - } - } + /* cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];or [codtabm(ij,TnsdVar[TvarsD[k]] */ + cov[2+nagesqr+k1]=Tresult[nres][resultmodel[nres][k1]]; + /* printf("hpxij Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,TnsdVar[TvarsD[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,TnsdVar[TvarsD[k]])],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,TnsdVar[TvarsD[k]])); */ + printf("hpxij Dummy combi=%d k1=%d Tvar[%d]=V%d cov[2+%d+%d]=%lf resultmodel[nres][%d]=%d nres/nresult=%d/%d \n",ij,k1,k1, Tvar[k1],nagesqr,k1,cov[2+nagesqr+k1],k1,resultmodel[nres][k1],nres,nresult); + }else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Single quantitative variables */ + /* resultmodel[nres][k1]=k3: k1th position in the model correspond to the k3 position in the resultline */ + cov[2+nagesqr+k1]=Tqresult[nres][resultmodel[nres][k1]]; + /* 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 k1=%d resultmodel[nres][%d]=%d,Tqresult[%d][%d]=%f\n",k1,k1,resultmodel[nres][k1],nres,resultmodel[nres][k1],Tqresult[nres][resultmodel[nres][k1]]); + }else if( Dummy[k1]==2 ){ /* For dummy with age product */ + /* Tvar[k1] Variable in the age product age*V1 is 1 */ + /* [Tinvresult[nres][V1] is its value in the resultline nres */ + cov[2+nagesqr+k1]=Tinvresult[nres][Tvar[k1]]; + printf("DhPxij Dummy with age k1=%d Tvar[%d]=%d Tinvresult[nres][%d]=%d,cov[2+%d+%d]=%.3f\n",k1,k1,Tvar[k1],Tinvresult[nres][Tvar[k1]],nagesqr,k1,cov[2+nagesqr+k1]); + /* cov[2+nagesqr+k1]=Tresult[nres][resultmodel[nres][k1]]; */ + /* for (k=1; k<=cptcovage;k++){ /\* For product with age V1+V1*age +V4 +age*V3 *\/ */ + /* 1+2 Tage[1]=2 TVar[2]=1 Dummy[2]=2, Tage[2]=4 TVar[4]=3 Dummy[4]=3 quant*/ + /* */ +/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ +/* k 1 2 3 4 5 6 7 8 9 */ +/*Tvar[k]= 5 4 3 6 5 2 7 1 1 */ +/*cptcovage=2 1 2 */ +/*Tage[k]= 5 8 */ + }else if( Dummy[k1]==3 ){ /* For quant with age product */ + cov[2+nagesqr+k1]=Tresult[nres][resultmodel[nres][k1]]; + printf("QhPxij Quant with age k1=%d resultmodel[nres][%d]=%d,Tqresult[%d][%d]=%f\n",k1,k1,resultmodel[nres][k1],nres,resultmodel[nres][k1],Tqresult[nres][resultmodel[nres][k1]]); + /* if(Dummy[Tage[k]]== 2){ /\* dummy with age *\/ */ + /* /\* if(Dummy[Tvar[Tage[k]]]== 2){ /\\* dummy with age *\\/ *\/ */ + /* /\* 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,TnsdVar[TvarsD[Tvar[Tage[k]]]])]*cov[2]; *\/ */ + /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,TnsdVar[TvarsD[Tvar[Tage[k]]]])]*cov[2]; */ + /* printf("hPxij Age combi=%d k=%d cptcovage=%d Tage[%d]=%d Tvar[Tage[%d]]=V%d nbcode[Tvar[Tage[k]]][codtabm(ij,TnsdVar[Tvar[Tage[k]]]])]=%d nres=%d\n",ij,k,cptcovage,k,Tage[k],k,Tvar[Tage[k]], nbcode[Tvar[Tage[k]]][codtabm(ij,TnsdVar[Tvar[Tage[k]]])],nres); */ + /* } else if(Dummy[Tage[k]]== 3){ /\* quantitative with age *\/ */ + /* 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]); */ + }else if(Typevar[k1]==2 ){ /* For product (not with age) */ +/* for (k=1; k<=cptcovprod;k++){ /\* For product without age *\/ */ +/* /\* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 *\/ */ +/* /\* k 1 2 3 4 5 6 7 8 9 *\/ */ +/* /\*Tvar[k]= 5 4 3 6 5 2 7 1 1 *\/ */ +/* /\*cptcovprod=1 1 2 *\/ */ +/* /\*Tprod[]= 4 7 *\/ */ +/* /\*Tvard[][1] 4 1 *\/ */ +/* /\*Tvard[][2] 3 2 *\/ */ + + /* printf("hPxij Prod ij=%d k=%d Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]=%d nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][1])]=%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2],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)]; */ + cov[2+nagesqr+k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]]; + printf("hPxij Prod ij=%d k1=%d cov[2+%d+%d]=%.5f Tvard[%d][1]=V%d * Tvard[%d][2]=V%d ; TinvDoQresult[nres][Tvardk[k1][1]]=%.4f * TinvDoQresult[nres][Tvardk[k1][1]]=%.4f\n",ij,k1,nagesqr,k1,cov[2+nagesqr+k1],k1,Tvard[k1][1], k1,Tvard[k1][2], TinvDoQresult[nres][Tvardk[k1][1]], TinvDoQresult[nres][Tvardk[k1][2]]); + /* if(Dummy[Tvardk[k1][1]]==0){ */ + /* if(Dummy[Tvardk[k1][2]]==0){ /\* Product of dummies *\/ */ + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */ + /* cov[2+nagesqr+k1]=Tinvresult[nres][Tvardk[k1][1]] * Tinvresult[nres][Tvardk[k1][2]]; */ + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,TnsdVar[Tvard[k][1]])] * nbcode[Tvard[k][2]][codtabm(ij,TnsdVar[Tvard[k][2]])]; */ + /* }else{ /\* Product of dummy by quantitative *\/ */ + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,TnsdVar[Tvard[k][1]])] * Tqresult[nres][k]; */ + /* cov[2+nagesqr+k1]=Tresult[nres][Tinvresult[nres][Tvardk[k1][1]]] * Tqresult[nres][Tinvresult[nres][Tvardk[k1][2]]]; */ + /* } */ + /* }else{ /\* Product of quantitative by...*\/ */ + /* if(Dummy[Tvard[k][2]]==0){ /\* quant by dummy *\/ */ + /* /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,TnsdVar[Tvard[k][2]])] * Tqinvresult[nres][Tvard[k][1]]; *\/ */ + /* cov[2+nagesqr+k1]=Tqresult[nres][Tinvresult[nres][Tvardk[k1][1]]] * Tresult[nres][Tinvresult[nres][Tvardk[k1][2]]] ; */ + /* }else{ /\* Product of two quant *\/ */ + /* /\* cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; *\/ */ + /* cov[2+nagesqr+k1]=Tqresult[nres][Tinvresult[nres][Tvardk[k1][1]]] * Tqresult[nres][Tinvresult[nres][Tvardk[k1][2]]] ; */ + /* } */ + /* }/\*end of products quantitative *\/ */ + }/*end of products */ + } /* End of loop on model equation */ /* 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 *\/ */ @@ -3550,7 +3619,7 @@ double ***hbxij(double ***po, int nhstep for (k=1; k<=nsd;k++){ /* For single dummy covariates only *//* cptcovn error */ /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */ /* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */ - cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];/* Bug valgrind */ + cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,TvarsD[k])];/* Bug valgrind */ /* printf("hbxij Dummy agexact=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agexact,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ } for (k=1; k<=nsq;k++) { /* For single varying covariates only */ @@ -3561,23 +3630,23 @@ double ***hbxij(double ***po, int nhstep for (k=1; k<=cptcovage;k++){ /* Should start at cptcovn+1 *//* For product with age */ /* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age error!!!*\/ */ if(Dummy[Tage[k]]== 2){ /* dummy with age */ - 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]; } else if(Dummy[Tage[k]]== 3){ /* quantitative with age */ cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; } /* printf("hBxij 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++){ /* Useless because included in cptcovn */ - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; - 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)]; + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,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,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][1])]; }else{ - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * 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]]; + if(Dummy[Tvard[k][2]]==0){ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])] * Tqinvresult[nres][Tvard[k][1]]; }else{ cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; } @@ -4790,6 +4859,7 @@ Title=%s
Datafile=%s Firstpass=%d La /* j=ncoveff; /\* Only fixed dummy covariates *\/ */ j=cptcoveff; /* Only dummy covariates of the model */ + /* j=cptcovn; /\* Only dummy covariates of the model *\/ */ if (cptcovn<1) {j=1;ncodemax[1]=1;} @@ -4797,7 +4867,7 @@ Title=%s
Datafile=%s Firstpass=%d La reference=low_education V1=0,V2=0 med_educ V1=1 V2=0, high_educ V1=0 V2=1 - Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff + Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcovn */ dateintsum=0; k2cpt=0; @@ -4809,7 +4879,7 @@ Title=%s
Datafile=%s Firstpass=%d La /* if a constant only model, one pass to compute frequency tables and to write it on ficresp */ /* Loop on nj=1 or 2 if dummy covariates j!=0 - * Loop on j1(1 to 2**cptcoveff) covariate combination + * Loop on j1(1 to 2**cptcovn) covariate combination * freq[s1][s2][iage] =0. * Loop on iind * ++freq[s1][s2][iage] weighted @@ -4834,12 +4904,12 @@ Title=%s
Datafile=%s Firstpass=%d La if(nj==1) j=0; /* First pass for the constant */ else{ - j=cptcoveff; /* Other passes for the covariate values */ + j=cptcovs; /* Other passes for the covariate values */ } first=1; for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on all covariates combination of the model, excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */ posproptt=0.; - /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); + /*printf("cptcovn=%d Tvaraff=%d", cptcovn,Tvaraff[1]); scanf("%d", i);*/ for (i=-5; i<=nlstate+ndeath; i++) for (s2=-5; s2<=nlstate+ndeath; s2++) @@ -4870,14 +4940,14 @@ Title=%s
Datafile=%s Firstpass=%d La bool=1; if(j !=0){ if(anyvaryingduminmodel==0){ /* If All fixed covariates */ - if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ - for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */ + if (cptcovn >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ + for (z1=1; z1<=cptcovn; z1++) { /* loops on covariates in the model */ /* if(Tvaraff[z1] ==-20){ */ /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ /* }else if(Tvaraff[z1] ==-10){ */ /* /\* sumnew+=coqvar[z1][iind]; *\/ */ - /* }else */ - if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ /* for combination j1 of covariates */ + /* }else */ /* TODO TODO codtabm(j1,z1) or codtabm(j1,Tvaraff[z1]]z1)*/ + if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]){ /* for combination j1 of covariates */ /* Tests if the value of the covariate z1 for this individual iind responded to combination j1 (V4=1 V3=0) */ bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */ /* 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", @@ -4895,15 +4965,15 @@ Title=%s
Datafile=%s Firstpass=%d La m=mw[mi][iind]; if(j!=0){ if(anyvaryingduminmodel==1){ /* Some are varying covariates */ - for (z1=1; z1<=cptcoveff; z1++) { + for (z1=1; z1<=cptcovn; z1++) { if( Fixed[Tmodelind[z1]]==1){ iv= Tvar[Tmodelind[z1]]-ncovcol-nqv; - if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality. If covariate's + if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]) /* iv=1 to ntv, right modality. If covariate's value is -1, we don't select. It differs from the constant and age model which counts them. */ bool=0; /* not selected */ }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */ - if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) { + if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]) { bool=0; } } @@ -4964,28 +5034,28 @@ Title=%s
Datafile=%s Firstpass=%d La /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ - if(cptcoveff==0 && nj==1) /* no covariate and first pass */ + if(cptcovn==0 && nj==1) /* no covariate and first pass */ pstamp(ficresp); - if (cptcoveff>0 && j!=0){ + if (cptcovn>0 && j!=0){ pstamp(ficresp); printf( "\n#********** Variable "); fprintf(ficresp, "\n#********** Variable "); fprintf(ficresphtm, "\n

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

********** Variable "); fprintf(ficlog, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++){ + for (z1=1; z1<=cptcovs; z1++){ if(!FixedV[Tvaraff[z1]]){ - printf( "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresp, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresphtm, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresphtmfr, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficlog, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + printf( "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); + fprintf(ficresp, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); + fprintf(ficresphtm, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); + fprintf(ficresphtmfr, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); + fprintf(ficlog, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); }else{ - printf( "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresphtm, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficresphtmfr, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); - fprintf(ficlog, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + printf( "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); + fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); + fprintf(ficresphtm, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); + fprintf(ficresphtmfr, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); + fprintf(ficlog, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); } } printf( "**********\n#"); @@ -5019,14 +5089,14 @@ Title=%s
Datafile=%s Firstpass=%d La /* } */ fprintf(ficresphtm,""); - if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */ + if((cptcovn==0 && nj==1)|| nj==2 ) /* no covariate and first pass */ fprintf(ficresp, " Age"); - if(nj==2) for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, " V%d=%d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + if(nj==2) for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " V%d=%d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]); for(i=1; i<=nlstate;i++) { - if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," Prev(%d) N(%d) N ",i,i); + if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp," Prev(%d) N(%d) N ",i,i); fprintf(ficresphtm, "",i,i); } - if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n"); + if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n"); fprintf(ficresphtm, "\n"); /* Header of frequency table by age */ @@ -5094,14 +5164,14 @@ Title=%s
Datafile=%s Firstpass=%d La } /* Writing ficresp */ - if(cptcoveff==0 && nj==1){ /* no covariate and first pass */ + if(cptcovn==0 && nj==1){ /* no covariate and first pass */ if( iage <= iagemax){ fprintf(ficresp," %d",iage); } }else if( nj==2){ if( iage <= iagemax){ fprintf(ficresp," %d",iage); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]); } } for(s1=1; s1 <=nlstate ; s1++){ @@ -5116,7 +5186,7 @@ Title=%s
Datafile=%s Firstpass=%d La } if( iage <= iagemax){ if(pos>=1.e-5){ - if(cptcoveff==0 && nj==1){ /* no covariate and first pass */ + if(cptcovn==0 && nj==1){ /* no covariate and first pass */ fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); }else if( nj==2){ fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); @@ -5125,7 +5195,7 @@ Title=%s
Datafile=%s Firstpass=%d La /*probs[iage][s1][j1]= pp[s1]/pos;*/ /*printf("\niage=%d s1=%d j1=%d %.5f %.0f %.0f %f",iage,s1,j1,pp[s1]/pos, pp[s1],pos,probs[iage][s1][j1]);*/ } else{ - if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta); + if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta); fprintf(ficresphtm,"",iage, prop[s1][iage],pospropta); } } @@ -5151,7 +5221,7 @@ Title=%s
Datafile=%s Firstpass=%d La } fprintf(ficresphtmfr,"\n "); fprintf(ficresphtm,"\n"); - if((cptcoveff==0 && nj==1)|| nj==2 ) { + if((cptcovn==0 && nj==1)|| nj==2 ) { if(iage <= iagemax) fprintf(ficresp,"\n"); } @@ -5427,10 +5497,10 @@ void prevalence(double ***probs, double for (z1=1; z1<=cptcoveff; z1++){ if( Fixed[Tmodelind[z1]]==1){ iv= Tvar[Tmodelind[z1]]-ncovcol-nqv; - if (cotvar[m][iv][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality */ + if (cotvar[m][iv][i]!= nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]) /* iv=1 to ntv, right modality */ bool=0; }else if( Fixed[Tmodelind[z1]]== 0) /* fixed */ - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) { + if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]) { bool=0; } } @@ -5876,7 +5946,7 @@ void concatwav(int wav[], int **dh, int } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */ /* ij--; */ /* cptcoveff=ij; /\*Number of total covariates*\/ */ - *cptcov=ij; /*Number of total real effective covariates: effective + *cptcov=ij; /* cptcov= Number of total real effective covariates: effective (used as cptcoveff in other functions) * because they can be excluded from the model and real * if in the model but excluded because missing values, but how to get k from ij?*/ for(j=ij+1; j<= cptcovt; j++){ @@ -5897,6 +5967,8 @@ void concatwav(int wav[], int **dh, int { /* Health expectancies, no variances */ + /* cij is the combination in the list of combination of dummy covariates */ + /* strstart is a string of time at start of computing */ int i, j, nhstepm, hstepm, h, nstepm; int nhstepma, nstepma; /* Decreasing with age */ double age, agelim, hf; @@ -5965,7 +6037,7 @@ void concatwav(int wav[], int **dh, int /* If stepm=6 months */ /* Computed by stepm unit matrices, product of hstepma matrices, stored in an array of nhstepma length: nhstepma=10, hstepm=4, stepm=6 months */ - + /* printf("HELLO evsij Entering hpxij age=%d cij=%d hstepm=%d x[1]=%f nres=%d\n",(int) age, cij, hstepm, x[1], nres); */ hpxij(p3mat,nhstepma,age,hstepm,x,nlstate,stepm,oldm, savm, cij, nres); hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ @@ -6279,7 +6351,7 @@ void concatwav(int wav[], int **dh, int 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,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(ij,Tvaraff[j])]); fprintf(ficresprobmorprev,"\n"); fprintf(ficresprobmorprev,"# Age cov=%-d",ij); @@ -6819,6 +6891,7 @@ void varprob(char optionfilefiname[], do int k2, l2, j1, z1; int k=0, l; int first=1, first1, first2; + int nres=0; /* New */ double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp; double **dnewm,**doldm; double *xp; @@ -6907,26 +6980,28 @@ To be simple, these graphs help to under if (cptcovn<1) {tj=1;ncodemax[1]=1;} j1=0; for(j1=1; j1<=tj;j1++){ /* For each valid combination of covariates or only once*/ + for(nres=1;nres <=1; nres++){ /* For each resultline */ + /* for(nres=1;nres <=nresult; nres++){ /\* For each resultline *\/ */ if (cptcovn>0) { fprintf(ficresprob, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[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]][codtabm(j1,z1)]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[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]][codtabm(j1,z1)]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[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]][codtabm(j1,z1)]); */ - for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtmcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtmcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]); fprintf(fichtmcov, "**********\n
"); fprintf(ficresprobcor, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,Tvaraff[z1])]); fprintf(ficresprobcor, "**********\n#"); if(invalidvarcomb[j1]){ fprintf(ficgp,"\n#Combination (%d) ignored because no cases \n",j1); @@ -6942,8 +7017,11 @@ To be simple, these graphs help to under cov[2]=age; if(nagesqr==1) cov[3]= age*age; - for (k=1; k<=cptcovn;k++) { - cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; + /* for (k=1; k<=cptcovn;k++) { */ + /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; */ + for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ + /* Here comes the value of the covariate 'j1' after renumbering k with single dummy covariates */ + cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(j1,TvarsD[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 @@ -6954,12 +7032,39 @@ To be simple, these graphs help to under /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */ /* ) p nbcode[Tvar[Tage[k]]][(1 & (ij-1) >> (k-1))+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]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; - for (k=1; k<=cptcovprod;k++) - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; - - + for (k=1; k<=cptcovage;k++){ /* For product with age */ + if(Dummy[Tage[k]]==2){ /* dummy with age */ + cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(j1,Tvar[Tage[k]])]*cov[2]; + /* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */ + } else if(Dummy[Tage[k]]==3){ /* quantitative with age */ + printf("Internal IMaCh error, don't know which value for quantitative covariate with age, Tage[k]%d, k=%d, Tvar[Tage[k]]=V%d, age=%d\n",Tage[k],k ,Tvar[Tage[k]], (int)cov[2]); + exit(1); + /* cov[2+nagesqr+Tage[k]]=meanq[k]/idq[k]*cov[2];/\* Using the mean of quantitative variable Tvar[Tage[k]] /\* Tqresult[nres][k]; *\/ */ + /* cov[++k1]=Tqresult[nres][k]; */ + } + /* cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */ + } + for (k=1; k<=cptcovprod;k++){/* For product without age */ + if(Dummy[Tvard[k][1]]==0){ + if(Dummy[Tvard[k][2]]==0){ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(j1,Tvard[k][2])]; + /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */ + }else{ /* Should we use the mean of the quantitative variables? */ + cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,Tvard[k][1])] * Tqresult[nres][k]; + /* cov[++k1]=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(j1,Tvard[k][2])] * Tqinvresult[nres][Tvard[k][1]]; + /* cov[++k1]=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]]; + /* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; */ + } + } + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; */ + } +/* For each age and combination of dummy covariates we slightly move the parameters of delti in order to get the gradient*/ for(theta=1; theta <=npar; theta++){ for(i=1; i<=npar; i++) xp[i] = x[i] + (i==theta ?delti[theta]:(double)0); @@ -7144,6 +7249,7 @@ To be simple, these graphs help to under } /* k12 */ } /*l1 */ }/* k1 */ + } /* loop on nres */ } /* loop on combination of covariates j1 */ free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); @@ -7294,19 +7400,23 @@ divided by h: hPij",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); /* Survival functions (period) in state j */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. %s_%d-%d-%d.svg
\ -", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); + fprintf(fichtm,"
\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. %s_%d-%d-%d.svg
", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); + fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_")); + fprintf(fichtm,"",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); } /* State specific survival functions (period) */ for(cpt=1; cpt<=nlstate;cpt++){ fprintf(fichtm,"
\n- Survival functions in state %d and in any other live state (total).\ And probability to be observed in various states (up to %d) being in state %d at different ages. \ - %s_%d-%d-%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); + %s_%d-%d-%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); + fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_")); + fprintf(fichtm,"",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); } /* Period (forward stable) prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. %s_%d-%d-%d.svg
\ -", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); + fprintf(fichtm,"
\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. %s_%d-%d-%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); + fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"P_"),subdirf2(optionfilefiname,"P_")); + fprintf(fichtm,"" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); } if(prevbcast==1){ /* Backward prevalence in each health state */ @@ -8367,41 +8477,49 @@ set ter svg size 640, 480\nunset log y\n /* for(j=3; j <=ncovmodel-nagesqr; j++) { */ for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */ /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */ - if(cptcovage >0){ /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */ - if(j==Tage[ij]) { /* Product by age To be looked at!!*//* Bug valgrind */ - if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */ - if(DummyV[j]==0){/* Bug valgrind */ - fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);; - }else{ /* quantitative */ - fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */ - /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */ + switch(Typevar[j]){ + case 1: + if(cptcovage >0){ /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */ + if(j==Tage[ij]) { /* Product by age To be looked at!!*//* Bug valgrind */ + if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */ + if(DummyV[j]==0){/* Bug valgrind */ + fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);; + }else{ /* quantitative */ + fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */ + /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */ + } + ij++; } - ij++; } - } - }else if(cptcovprod >0){ - if(j==Tprod[ijp]) { /* */ - /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */ - if(ijp <=cptcovprod) { /* Product */ - if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */ - if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */ - /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */ - fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]); - }else{ /* Vn is dummy and Vm is quanti */ - /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */ - fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); - } - }else{ /* Vn*Vm Vn is quanti */ - if(DummyV[Tvard[ijp][2]]==0){ - fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]); - }else{ /* Both quanti */ - fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); + } + break; + case 2: + if(cptcovprod >0){ + if(j==Tprod[ijp]) { /* */ + /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */ + if(ijp <=cptcovprod) { /* Product */ + if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */ + if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */ + /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */ + fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]); + }else{ /* Vn is dummy and Vm is quanti */ + /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */ + fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); + } + }else{ /* Vn*Vm Vn is quanti */ + if(DummyV[Tvard[ijp][2]]==0){ + fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]); + }else{ /* Both quanti */ + fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); + } } + ijp++; } - ijp++; - } - } /* end Tprod */ - } else{ /* simple covariate */ + } /* end Tprod */ + } + break; + case 0: + /* simple covariate */ /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(k1,j)]); /\* Valgrind bug nbcode *\/ */ if(Dummy[j]==0){ fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]); /* */ @@ -8409,12 +8527,17 @@ set ter svg size 640, 480\nunset log y\n fprintf(ficgp,"+p%d*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* */ /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */ } - } /* end simple */ + /* end simple */ + break; + default: + break; + } /* end switch */ } /* end j */ - }else{ - i=i-ncovmodel; - if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */ - fprintf(ficgp," (1."); + }else{ /* k=k2 */ + if(ng !=1 ){ /* For logit formula of log p11 is more difficult to get */ + fprintf(ficgp," (1.");i=i-ncovmodel; + }else + i=i-ncovmodel; } if(ng != 1){ @@ -8427,17 +8550,78 @@ set ter svg size 640, 480\nunset log y\n fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(cpt-1)*ncovmodel,k3+(cpt-1)*ncovmodel+1,k3+(cpt-1)*ncovmodel+1+nagesqr); ij=1; - for(j=3; j <=ncovmodel-nagesqr; j++){ - if(cptcovage >0){ - if((j-2)==Tage[ij]) { /* Bug valgrind */ - if(ij <=cptcovage) { /* Bug valgrind */ - fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]); - /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */ - ij++; - } - } - }else - fprintf(ficgp,"+p%d*%d",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);/* Valgrind bug nbcode */ + ijp=1; + /* for(j=3; j <=ncovmodel-nagesqr; j++){ */ + for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */ + switch(Typevar[j]){ + case 1: + if(cptcovage >0){ + if(j==Tage[ij]) { /* Bug valgrind */ + if(ij <=cptcovage) { /* Bug valgrind */ + if(DummyV[j]==0){/* Bug valgrind */ + /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]); */ + /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,nbcode[Tvar[j]][codtabm(k1,j)]); */ + fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvar[j]]); + /* fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);; */ + /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */ + }else{ /* quantitative */ + /* fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /\* Tqinvresult in decoderesult *\/ */ + fprintf(ficgp,"+p%d*%f*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */ + /* fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /\* Tqinvresult in decoderesult *\/ */ + /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */ + } + ij++; + } + } + } + break; + case 2: + if(cptcovprod >0){ + if(j==Tprod[ijp]) { /* */ + /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */ + if(ijp <=cptcovprod) { /* Product */ + if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */ + if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */ + /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */ + fprintf(ficgp,"+p%d*%d*%d",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]); + /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]); */ + }else{ /* Vn is dummy and Vm is quanti */ + /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */ + fprintf(ficgp,"+p%d*%d*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); + /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); */ + } + }else{ /* Vn*Vm Vn is quanti */ + if(DummyV[Tvard[ijp][2]]==0){ + fprintf(ficgp,"+p%d*%d*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]); + /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]); */ + }else{ /* Both quanti */ + fprintf(ficgp,"+p%d*%f*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); + /* fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); */ + } + } + ijp++; + } + } /* end Tprod */ + } /* end if */ + break; + case 0: + /* simple covariate */ + /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(k1,j)]); /\* Valgrind bug nbcode *\/ */ + if(Dummy[j]==0){ + /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]); /\* *\/ */ + fprintf(ficgp,"+p%d*%d",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvar[j]]); /* */ + /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]); /\* *\/ */ + }else{ /* quantitative */ + fprintf(ficgp,"+p%d*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tqinvresult[nres][Tvar[j]]); /* */ + /* fprintf(ficgp,"+p%d*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /\* *\/ */ + /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */ + } + /* end simple */ + /* fprintf(ficgp,"+p%d*%d",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);/\* Valgrind bug nbcode *\/ */ + break; + default: + break; + } /* end switch */ } fprintf(ficgp,")"); } @@ -8446,7 +8630,7 @@ set ter svg size 640, 480\nunset log y\n fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"p%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k); else /* ng= 3 */ fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"i%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k); - }else{ /* end ng <> 1 */ + }else{ /* end ng <> 1 */ if( k !=k2) /* logit p11 is hard to draw */ fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"logit(p%d%d)\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k); } @@ -8791,7 +8975,7 @@ void prevforecast(char fileres[], double } fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#"); for(j=1;j<=cptcoveff;j++) { - fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); } for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); @@ -8821,7 +9005,7 @@ void prevforecast(char fileres[], double } fprintf(ficresf,"\n"); for(j=1;j<=cptcoveff;j++) - fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); fprintf(ficresf,"%.f %.f ",anprojd+yearp,agec+h*hstepm/YEARM*stepm); for(j=1; j<=nlstate+ndeath;j++) { @@ -8931,7 +9115,7 @@ void prevforecast(char fileres[], double } fprintf(ficresfb,"\n#****** hbijx=probability over h years, hb.jx is weighted by observed prev \n#"); for(j=1;j<=cptcoveff;j++) { - fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); } for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); @@ -8967,7 +9151,7 @@ void prevforecast(char fileres[], double } fprintf(ficresfb,"\n"); for(j=1;j<=cptcoveff;j++) - fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); fprintf(ficresfb,"%.f %.f ",anbackd+yearp,agec-h*hstepm/YEARM*stepm); for(i=1; i<=nlstate+ndeath;i++) { ppij=0.;ppi=0.; @@ -9034,9 +9218,9 @@ void prevforecast(char fileres[], double printf("\n#****** "); fprintf(ficlog,"\n#****** "); for(j=1;j<=cptcoveff;j++) { - fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); + fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); + printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); } for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); @@ -9091,9 +9275,9 @@ void prevforecast(char fileres[], double printf("\n#****** "); fprintf(ficlog,"\n#****** "); for(j=1;j<=cptcoveff;j++) { - fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); + fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); + printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); } for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); @@ -9560,6 +9744,10 @@ int readdata(char datafile[], int firsto DummyV=ivector(1,NCOVMAX); /* 1 to 3 */ FixedV=ivector(1,NCOVMAX); /* 1 to 3 */ + for(v=1;v V4*2**(0) + V3*2**(1) + V1*2**(2) = 5 + (1offset) = 6*/ - /* result line V4=1 V5=24.1 V3=1 V2=8 V1=0 */ + /* nres=1st result line: V4=1 V5=25.1 V3=0 V2=8 V1=1 */ + /* should correspond to the combination 6 of dummy: V4=1, V3=0, V1=1 => V4*2**(0) + V3*2**(1) + V1*2**(2) = 1*1 + 0*2 + 1*4 = 5 + (1offset) = 6*/ + /* nres=2nd result line: V4=1 V5=24.1 V3=1 V2=8 V1=0 */ /* should give a combination of dummy V4=1, V3=1, V1=0 => V4*2**(0) + V3*2**(1) + V1*2**(2) = 3 + (1offset) = 4*/ /* 1 0 0 0 */ /* 2 1 0 0 */ /* 3 0 1 0 */ - /* 4 1 1 0 */ /* V4=1, V3=1, V1=0 */ + /* 4 1 1 0 */ /* V4=1, V3=1, V1=0 (nres=2)*/ /* 5 0 0 1 */ - /* 6 1 0 1 */ /* V4=1, V3=0, V1=1 */ + /* 6 1 0 1 */ /* V4=1, V3=0, V1=1 (nres=1)*/ /* 7 0 1 1 */ /* 8 1 1 1 */ /* V(Tvresult)=Tresult V4=1 V3=0 V1=1 Tresult[nres=1][2]=0 */ /* V(Tvqresult)=Tqresult V5=25.1 V2=8 Tqresult[nres=1][1]=25.1 */ /* V5*age V5 known which value for nres? */ /* Tqinvresult[2]=8 Tqinvresult[1]=25.1 */ - for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* loop on model line */ - if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Single dummy */ - k3= resultmodel[k1]; /* resultmodel[2(V4)] = 1=k3 */ - k2=(int)Tvarsel[k3]; /* Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */ - k+=Tvalsel[k3]*pow(2,k4); /* Tvalsel[1]=1 */ - Tresult[nres][k4+1]=Tvalsel[k3];/* Tresult[nres][1]=1(V4=1) Tresult[nres][2]=0(V3=0) */ + for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* loop k1 on position in the model line (excluding product) */ + /* k counting number of combination of single dummies in the equation model */ + /* k4 counting single dummies in the equation model */ + /* k4q counting single quantitatives in the equation model */ + if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Dummy and Single */ + /* k4+1= position in the resultline V(Tvarsel)=Tvalsel=Tresult[nres][pos](value); V(Tvresult[nres][pos] (variable): V(variable)=value) */ + /* modelresult[k3]=k1: k3th position in the result line correspond to the k1 position in the model line */ + /* Value in the (current nres) resultline of the variable at the k1th position in the model equation resultmodel[nres][k1]= k3 */ + /* resultmodel[nres][k1]=k3: k1th position in the model correspond to the k3 position in the resultline */ + /* k3 is the position in the nres result line of the k1th variable of the model equation */ + /* Tvarsel[k3]: Name of the variable at the k3th position in the result line. */ + /* Tvalsel[k3]: Value of the variable at the k3th position in the result line. */ + /* Tresult[nres][result_position]= value of the dummy variable at the result_position in the nres resultline */ + /* Tvresult[nres][result_position]= id of the dummy variable at the result_position in the nres resultline */ + /* Tinvresult[nres][Name of a dummy variable]= value of the variable in the result line */ + /* TinvDoQresult[nres][Name of a Dummy or Q variable]= value of the variable in the result line */ + k3= resultmodel[nres][k1]; /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/ + k2=(int)Tvarsel[k3]; /* nres=1 k1=2=>k3=1 Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 (V4); k1=3=>k3=2 Tvarsel[2]=3 (V3)*/ + k+=Tvalsel[k3]*pow(2,k4); /* nres=1 k1=2 Tvalsel[1]=1 (V4=1); k1=3 k3=2 Tvalsel[2]=0 (V3=0) */ + TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */ + Tresult[nres][k4+1]=Tvalsel[k3];/* Tresult[nres=2][1]=1(V4=1) Tresult[nres=2][2]=0(V3=0) */ Tvresult[nres][k4+1]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */ Tinvresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */ printf("Decoderesult Dummy k=%d, V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k2, k3, (int)Tvalsel[k3], k4); k4++;; - } else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Single quantitative */ - k3q= resultmodel[k1]; /* resultmodel[1(V5)] = 25.1=k3q */ + }else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Quantitative and single */ + /* Tqresult[nres][result_position]= value of the variable at the result_position in the nres resultline */ + /* Tqvresult[nres][result_position]= id of the variable at the result_position in the nres resultline */ + /* Tqinvresult[nres][Name of a quantitative variable]= value of the variable in the result line */ + k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 25.1=k3q */ k2q=(int)Tvarsel[k3q]; /* Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */ Tqresult[nres][k4q+1]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */ Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */ Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */ + TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */ printf("Decoderesult Quantitative nres=%d, V(k2q=V%d)= Tvalsel[%d]=%d, Tvarsel[%d]=%f\n",nres, k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]); k4q++;; + }else if( Dummy[k1]==2 ){ /* For dummy with age product */ + /* Tvar[k1]; */ /* Age variable */ + k3= resultmodel[nres][Tvar[k1]]; /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/ + k2=(int)Tvarsel[k3]; /* nres=1 k1=2=>k3=1 Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 (V4); k1=3=>k3=2 Tvarsel[2]=3 (V3)*/ + TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */ + printf("Decoderesult Dummy with age k=%d, k1=%d Tvar[%d]=V%d k2=Tvarsel[%d]=%d Tvalsel[%d]=%d\n",k, k1, k1, Tvar[k1], k3,(int)Tvarsel[k3], k3, (int)Tvalsel[k3]); + }else if( Dummy[k1]==3 ){ /* For quant with age product */ + k3q= resultmodel[nres][Tvar[k1]]; /* resultmodel[1(V5)] = 25.1=k3q */ + k2q=(int)Tvarsel[k3q]; /* Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */ + TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */ + printf("Decoderesult Quantitative with age nres=%d, k1=%d, Tvar[%d]=V%d V(k2q=V%d)= Tvalsel[%d]=%d, Tvarsel[%d]=%f\n",nres, k1, k1, Tvar[k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]); + }else if(Typevar[k1]==2 ){ /* For product quant or dummy (not with age) */ + printf("Decoderesult Quantitative or Dummy (not with age) nres=%d k1=%d Tvar[%d]=%d \n",nres, k1, k1, Tvar[k1]); + }else{ + printf("Error Decodemodel probably a product Dummy[%d]==%d && Typevar[%d]==%d\n", k1, Dummy[k1], k1, Typevar[k1]); + fprintf(ficlog,"Error Decodemodel probably a product Dummy[%d]==%d && Typevar[%d]==%d\n", k1, Dummy[k1], k1, Typevar[k1]); } } @@ -10209,7 +10432,9 @@ int decodemodel( char model[], int lasto Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */ Tposprod[k]=k1; /* Tposprod[3]=1, Tposprod[2]=5 */ Tvard[k1][1] =atoi(strc); /* m 1 for V1*/ + Tvardk[k][1] =atoi(strc); /* m 1 for V1*/ Tvard[k1][2] =atoi(stre); /* n 4 for V4*/ + Tvardk[k][2] =atoi(stre); /* n 4 for V4*/ k2=k2+2; /* k2 is initialize to -1, We want to store the n and m in Vn*Vm at the end of Tvar */ /* Tvar[cptcovt+k2]=Tvard[k1][1]; /\* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) *\/ */ /* Tvar[cptcovt+k2+1]=Tvard[k1][2]; /\* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) *\/ */ @@ -10281,6 +10506,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( modell[k].maintype= FTYPE; TvarsD[nsd]=Tvar[k]; TvarsDind[nsd]=k; + TnsdVar[Tvar[k]]=nsd; TvarF[ncovf]=Tvar[k]; TvarFind[ncovf]=k; TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ @@ -10292,6 +10518,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( ncovf++; modell[k].maintype= FTYPE; TvarF[ncovf]=Tvar[k]; + /* TnsdVar[Tvar[k]]=nsd; */ /* To be done */ TvarFind[ncovf]=k; TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ @@ -10318,6 +10545,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( nsd++; TvarsD[nsd]=Tvar[k]; TvarsDind[nsd]=k; + TnsdVar[Tvar[k]]=nsd; /* To be verified */ ncovv++; /* Only simple time varying variables */ TvarV[ncovv]=Tvar[k]; TvarVind[ncovv]=k; /* TvarVind[2]=2 TvarVind[3]=3 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */ @@ -10903,9 +11131,9 @@ int prevalence_limit(double *p, double * printf("#******"); fprintf(ficlog,"#******"); for(j=1;j<=cptcoveff ;j++) {/* all covariates */ - fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); /* Here problem for varying dummy*/ - printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); /* Here problem for varying dummy*/ + printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); + fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); } for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); @@ -10924,7 +11152,7 @@ int prevalence_limit(double *p, double * fprintf(ficrespl,"#Age "); for(j=1;j<=cptcoveff;j++) { - fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); } for(i=1; i<=nlstate;i++) fprintf(ficrespl," %d-%d ",i,i); fprintf(ficrespl,"Total Years_to_converge\n"); @@ -10934,7 +11162,7 @@ int prevalence_limit(double *p, double * prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k, nres); fprintf(ficrespl,"%.0f ",age ); for(j=1;j<=cptcoveff;j++) - fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); tot=0.; for(i=1; i<=nlstate;i++){ tot += prlim[i][i]; @@ -10994,9 +11222,9 @@ int back_prevalence_limit(double *p, dou printf("#******"); fprintf(ficlog,"#******"); for(j=1;j<=cptcoveff ;j++) {/* all covariates */ - fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); + printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); + fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); } for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); @@ -11015,7 +11243,7 @@ int back_prevalence_limit(double *p, dou fprintf(ficresplb,"#Age "); for(j=1;j<=cptcoveff;j++) { - fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); } for(i=1; i<=nlstate;i++) fprintf(ficresplb," %d-%d ",i,i); fprintf(ficresplb,"Total Years_to_converge\n"); @@ -11039,7 +11267,7 @@ int back_prevalence_limit(double *p, dou } fprintf(ficresplb,"%.0f ",age ); for(j=1;j<=cptcoveff;j++) - fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); tot=0.; for(i=1; i<=nlstate;i++){ tot += bprlim[i][i]; @@ -11097,7 +11325,7 @@ int hPijx(double *p, int bage, int fage) continue; fprintf(ficrespij,"\n#****** "); for(j=1;j<=cptcoveff;j++) - fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); fprintf(ficrespij," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); @@ -11176,7 +11404,7 @@ int hPijx(double *p, int bage, int fage) continue; fprintf(ficrespijb,"\n#****** "); for(j=1;j<=cptcoveff;j++) - fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ fprintf(ficrespijb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); } @@ -11917,6 +12145,7 @@ Please run with mle=-1 to get a correct Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */ TvarsDind=ivector(1,NCOVMAX); /* */ + TnsdVar=ivector(1,NCOVMAX); /* */ TvarsD=ivector(1,NCOVMAX); /* */ TvarsQind=ivector(1,NCOVMAX); /* */ TvarsQ=ivector(1,NCOVMAX); /* */ @@ -11959,10 +12188,13 @@ Please run with mle=-1 to get a correct Tvard=imatrix(1,NCOVMAX,1,2); /* n=Tvard[k1][1] and m=Tvard[k1][2] gives the couple n,m of the k1 th product Vn*Vm * For V3*V2 (in V2+V1+V1*V4+age*V3+V3*V2), V3*V2 position is 2nd. * Tvard[k1=2][1]=3 (V3) Tvard[k1=2][2]=2(V2) */ + Tvardk=imatrix(1,NCOVMAX,1,2); Tage=ivector(1,NCOVMAX); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age 4 covariates (3 plus signs) Tage[1=V3*age]= 4; Tage[2=age*V4] = 3 - */ + */ + for(i=1;i> (k-1)) + 1; * #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 - * h\k 1 2 3 4 - *______________________________ - * 1 i=1 1 i=1 1 i=1 1 i=1 1 - * 2 2 1 1 1 - * 3 i=2 1 2 1 1 - * 4 2 2 1 1 - * 5 i=3 1 i=2 1 2 1 - * 6 2 1 2 1 - * 7 i=4 1 2 2 1 - * 8 2 2 2 1 - * 9 i=5 1 i=3 1 i=2 1 2 - * 10 2 1 1 2 - * 11 i=6 1 2 1 2 - * 12 2 2 1 2 - * 13 i=7 1 i=4 1 2 2 - * 14 2 1 2 2 - * 15 i=8 1 2 2 2 - * 16 2 2 2 2 - */ + * h\k 1 2 3 4 * h-1\k-1 4 3 2 1 + *______________________________ *______________________ + * 1 i=1 1 i=1 1 i=1 1 i=1 1 * 0 0 0 0 0 + * 2 2 1 1 1 * 1 0 0 0 1 + * 3 i=2 1 2 1 1 * 2 0 0 1 0 + * 4 2 2 1 1 * 3 0 0 1 1 + * 5 i=3 1 i=2 1 2 1 * 4 0 1 0 0 + * 6 2 1 2 1 * 5 0 1 0 1 + * 7 i=4 1 2 2 1 * 6 0 1 1 0 + * 8 2 2 2 1 * 7 0 1 1 1 + * 9 i=5 1 i=3 1 i=2 1 2 * 8 1 0 0 0 + * 10 2 1 1 2 * 9 1 0 0 1 + * 11 i=6 1 2 1 2 * 10 1 0 1 0 + * 12 2 2 1 2 * 11 1 0 1 1 + * 13 i=7 1 i=4 1 2 2 * 12 1 1 0 0 + * 14 2 1 2 2 * 13 1 1 0 1 + * 15 i=8 1 2 2 2 * 14 1 1 1 0 + * 16 2 2 2 2 * 15 1 1 1 1 + */ /* How to do the opposite? From combination h (=1 to 2**k) how to get the value on the covariates? */ /* from h=5 and m, we get then number of covariates k=log(m)/log(2)=4 * and the value of each covariate? @@ -13160,8 +13392,8 @@ Please run with mle=-1 to get a correct fprintf(ficreseij,"\n#****** "); printf("\n#****** "); for(j=1;j<=cptcoveff;j++) { - fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); + printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); } for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); @@ -13172,6 +13404,7 @@ Please run with mle=-1 to get a correct eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); oldm=oldms;savm=savms; + /* printf("HELLO Entering evsij bage=%d fage=%d k=%d estepm=%d nres=%d\n",(int) bage, (int)fage, k, estepm, nres); */ evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart, nres); free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); @@ -13229,9 +13462,9 @@ Please run with mle=-1 to get a correct fprintf(ficrest,"\n# model %s \n#****** Result for:", model); fprintf(ficlog,"\n# model %s \n#****** Result for:", model); for(j=1;j<=cptcoveff;j++){ - printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); + fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); + fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); } for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); @@ -13245,8 +13478,8 @@ Please run with mle=-1 to get a correct fprintf(ficresstdeij,"\n#****** "); fprintf(ficrescveij,"\n#****** "); for(j=1;j<=cptcoveff;j++) { - fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); - fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); + fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); } for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ fprintf(ficresstdeij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); @@ -13258,7 +13491,7 @@ Please run with mle=-1 to get a correct fprintf(ficresvij,"\n#****** "); /* pstamp(ficresvij); */ for(j=1;j<=cptcoveff;j++) - fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); } @@ -13348,6 +13581,7 @@ Please run with mle=-1 to get a correct free_vector(weight,firstobs,lastobs); + free_imatrix(Tvardk,1,NCOVMAX,1,2); free_imatrix(Tvard,1,NCOVMAX,1,2); free_imatrix(s,1,maxwav+1,firstobs,lastobs); free_matrix(anint,1,maxwav,firstobs,lastobs); @@ -13395,6 +13629,7 @@ Please run with mle=-1 to get a correct free_ivector(TvarsQ,1,NCOVMAX); free_ivector(TvarsQind,1,NCOVMAX); free_ivector(TvarsD,1,NCOVMAX); + free_ivector(TnsdVar,1,NCOVMAX); free_ivector(TvarsDind,1,NCOVMAX); free_ivector(TvarFD,1,NCOVMAX); free_ivector(TvarFDind,1,NCOVMAX);
AgePrev(%d)N(%d)N%dNaNq%.0f%.0f