Summary: r30

* imach.c (Module): Error cptcovn instead of nsd in bmij (was
coredumped, revealed by Feiuno, thank you.

diff --git a/src/imach.c b/src/imach.c
index 6e7a58d6bd0af69205389f402e0d93a3bdc96e2e..cb9e39cff7f5d4eff7c2ced7c1d1d69884c9f65b 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -1,6 +1,9 @@
 /* $Id$
   $State$
   $Log$
+  Revision 1.324  2022/07/23 17:44:26  brouard
+  *** empty log message ***
+
   Revision 1.323  2022/07/22 12:30:08  brouard
   *  imach.c (Module): Output of Wald test in the htm file and not only in the log.
 
@@ -1048,7 +1051,7 @@ Important routines
 - Tricode which tests the modality of dummy variables (in order to warn with wrong or empty modalities)
   and returns the number of efficient covariates cptcoveff and modalities nbcode[Tvar[k]][1]= 0 and nbcode[Tvar[k]][2]= 1 usually.
 - printinghtml which outputs results like life expectancy in and from a state for a combination of modalities of dummy variables
-  o There are 2*cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, éliminating 1 1 if
+  o There are 2**cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, éliminating 1 1 if
     race White (0 0), Black vs White (1 0), Hispanic (0 1) and 1 1 being meaningless.
 
 
@@ -1185,7 +1188,7 @@ typedef struct {
 #define NINTERVMAX 8
 #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */
 #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
-#define NCOVMAX 30  /**< Maximum number of covariates, including generated covariates V1*V2 */
+#define NCOVMAX 30  /**< Maximum number of covariates used in the model, including generated covariates V1*V2 or V1*age */
 #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
 /*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/
 #define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 
@@ -3135,7 +3138,7 @@ double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
     ps[i][i]=1./(s1+1.);
     /* Computing other pijs */
     for(j=1; j<i; j++)
-      ps[i][j]= exp(ps[i][j])*ps[i][i];
+      ps[i][j]= exp(ps[i][j])*ps[i][i];/* Bug valgrind */
     for(j=i+1; j<=nlstate+ndeath; j++)
       ps[i][j]= exp(ps[i][j])*ps[i][i];
     /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
@@ -3193,7 +3196,7 @@ double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
   /* dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
 
   /* P_x */
-  pmmij=pmij(pmmij,cov,ncovmodel,x,nlstate); /*This is forward probability from agefin to agefin + stepm */
+  pmmij=pmij(pmmij,cov,ncovmodel,x,nlstate); /*This is forward probability from agefin to agefin + stepm *//* Bug valgrind */
   /* outputs pmmij which is a stochastic matrix in row */
 
   /* Diag(w_x) */
@@ -3538,10 +3541,10 @@ double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, do
       cov[2]=agexact;
       if(nagesqr==1)
        cov[3]= agexact*agexact;
-      for (k=1; k<=cptcovn;k++){
+      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)];
+       cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,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 */
@@ -3560,6 +3563,19 @@ double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, do
       }
       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)];
+         }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("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
       /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
@@ -3569,7 +3585,7 @@ double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, do
       /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */
       /*                                                1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */
       out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\
-                  1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
+                  1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);/* Bug valgrind */
       /* if((int)age == 70){ */
       /*       printf(" Backward hbxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */
       /*       for(i=1; i<=nlstate+ndeath; i++) { */
@@ -8346,9 +8362,9 @@ set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar)
            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!!*/
+               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){
+                   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 */
@@ -11177,7 +11193,7 @@ int hPijx(double *p, int bage, int fage){
 
        /* oldm=oldms;savm=savms; */
        /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);   */
-       hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k, nres);
+       hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k, nres);/* Bug valgrind */
        /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
        fprintf(ficrespijb,"# Cov Agex agex-h hbijx with i,j=");
        for(i=1; i<=nlstate;i++)
@@ -11190,7 +11206,7 @@ int hPijx(double *p, int bage, int fage){
          /* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */
          for(i=1; i<=nlstate;i++)
            for(j=1; j<=nlstate+ndeath;j++)
-             fprintf(ficrespijb," %.5f", p3mat[i][j][h]);
+             fprintf(ficrespijb," %.5f", p3mat[i][j][h]);/* Bug valgrind */
          fprintf(ficrespijb,"\n");
        }
        free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
@@ -11875,7 +11891,8 @@ Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numline
   }
   mint=matrix(1,maxwav,firstobs,lastobs);
   anint=matrix(1,maxwav,firstobs,lastobs);
-  s=imatrix(1,maxwav+1,firstobs,lastobs); /* s[i][j] health state for wave i and individual j */ 
+  s=imatrix(1,maxwav+1,firstobs,lastobs); /* s[i][j] health state for wave i and individual j */
+  printf("BUG ncovmodel=%d NCOVMAX=%d 2**ncovmodel=%f BUG\n",ncovmodel,NCOVMAX,pow(2,ncovmodel));
   tab=ivector(1,NCOVMAX);
   ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */
   ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */