[BACK]Return to imach.c CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/imach.c between versions 1.319 and 1.329

version 1.319, 2022/06/02 04:45:11 version 1.329, 2022/08/03 17:29:54
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     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
   
     * imach.c (Module): Error cptcovn instead of nsd in bmij (was
     coredumped, revealed by Feiuno, thank you.
   
     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.
   
     Revision 1.322  2022/07/22 12:27:48  brouard
     *  imach.c (Module): Output of Wald test in the htm file and not only in the log.
   
     Revision 1.321  2022/07/22 12:04:24  brouard
     Summary: r28
   
     *  imach.c (Module): Output of Wald test in the htm file and not only in the log.
   
     Revision 1.320  2022/06/02 05:10:11  brouard
     *** empty log message ***
   
   Revision 1.319  2022/06/02 04:45:11  brouard    Revision 1.319  2022/06/02 04:45:11  brouard
   * imach.c (Module): Adding the Wald tests from the log to the main    * imach.c (Module): Adding the Wald tests from the log to the main
   htm for better display of the maximum likelihood estimators.    htm for better display of the maximum likelihood estimators.
Line 855 Line 890
   
   The same imach parameter file can be used but the option for mle should be -3.    The same imach parameter file can be used but the option for mle should be -3.
   
   Agnès, who wrote this part of the code, tried to keep most of the    Agnès, who wrote this part of the code, tried to keep most of the
   former routines in order to include the new code within the former code.    former routines in order to include the new code within the former code.
   
   The output is very simple: only an estimate of the intercept and of    The output is very simple: only an estimate of the intercept and of
Line 1034  Important routines Line 1069  Important routines
 - Tricode which tests the modality of dummy variables (in order to warn with wrong or empty modalities)  - 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.    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  - 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.      race White (0 0), Black vs White (1 0), Hispanic (0 1) and 1 1 being meaningless.
   
   
       
   Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).    Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).
            Institut national d'études démographiques, Paris.             Institut national d'études démographiques, Paris.
   This software have been partly granted by Euro-REVES, a concerted action    This software have been partly granted by Euro-REVES, a concerted action
   from the European Union.    from the European Union.
   It is copyrighted identically to a GNU software product, ie programme and    It is copyrighted identically to a GNU software product, ie programme and
Line 1160  typedef struct { Line 1195  typedef struct {
   
 #define GNUPLOTPROGRAM "gnuplot"  #define GNUPLOTPROGRAM "gnuplot"
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/  /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
 #define FILENAMELENGTH 132  #define FILENAMELENGTH 256
   
 #define GLOCK_ERROR_NOPATH              -1      /* empty path */  #define GLOCK_ERROR_NOPATH              -1      /* empty path */
 #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */  #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */
Line 1171  typedef struct { Line 1206  typedef struct {
 #define NINTERVMAX 8  #define NINTERVMAX 8
 #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */  #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */
 #define NDEATHMAX 8 /**< Maximum number of dead 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 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<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/
 #define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1   #define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 
Line 1199  typedef struct { Line 1234  typedef struct {
 /* $State$ */  /* $State$ */
 #include "version.h"  #include "version.h"
 char version[]=__IMACH_VERSION__;  char version[]=__IMACH_VERSION__;
 char copyright[]="May 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 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$ $Date$";   char fullversion[]="$Revision$ $Date$"; 
 char strstart[80];  char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];  char optionfilext[10], optionfilefiname[FILENAMELENGTH];
Line 1422  int **nbcode, *Tvar; /**< model=V2 => Tv Line 1457  int **nbcode, *Tvar; /**< model=V2 => Tv
 /* Tage[cptcovage]=k            5               8      */ /* Position in the model of ith cov*age */  /* 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 */  /* Tvard[1][1]@4={4,3,1,2}    V4*V3 V1*V2              */ /* Position in model of the ith prod without age */
 /* 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 */  /* 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)  */  /* TvarFind;  TvarFind[1]=6,  TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod)  */
 /* Type                    */  /* Type                    */
 /* V         1  2  3  4  5 */  /* V         1  2  3  4  5 */
 /*           F  F  V  V  V */  /*           F  F  V  V  V */
Line 2402  void powell(double p[], double **xi, int Line 2437  void powell(double p[], double **xi, int
   for (j=1;j<=n;j++) pt[j]=p[j];     for (j=1;j<=n;j++) pt[j]=p[j]; 
   rcurr_time = time(NULL);      rcurr_time = time(NULL);  
   for (*iter=1;;++(*iter)) {     for (*iter=1;;++(*iter)) { 
     fp=(*fret); /* From former iteration or initial value */  
     ibig=0;       ibig=0; 
     del=0.0;       del=0.0; 
     rlast_time=rcurr_time;      rlast_time=rcurr_time;
     /* (void) gettimeofday(&curr_time,&tzp); */      /* (void) gettimeofday(&curr_time,&tzp); */
     rcurr_time = time(NULL);        rcurr_time = time(NULL);  
     curr_time = *localtime(&rcurr_time);      curr_time = *localtime(&rcurr_time);
     printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);      printf("\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);
     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(ficlog,"\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*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); */  /*     fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */
       fp=(*fret); /* From former iteration or initial value */
     for (i=1;i<=n;i++) {      for (i=1;i<=n;i++) {
       fprintf(ficrespow," %.12lf", p[i]);        fprintf(ficrespow," %.12lf", p[i]);
     }      }
Line 2792  void powell(double p[], double **xi, int Line 2827  void powell(double p[], double **xi, int
     }      }
     for (k=1; k<=cptcovprod;k++){ /* For product without age */      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]); */        /* 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][1]]==0){
         if(Dummy[Tvard[k][2]==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,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
           /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */            /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
         }else{          }else{
Line 2801  void powell(double p[], double **xi, int Line 2836  void powell(double p[], double **xi, int
           /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; */            /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; */
         }          }
       }else{        }else{
         if(Dummy[Tvard[k][2]==0]){          if(Dummy[Tvard[k][2]]==0){
           cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
           /* cov[++k1]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; */            /* cov[++k1]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; */
         }else{          }else{
Line 2971  void powell(double p[], double **xi, int Line 3006  void powell(double p[], double **xi, int
     }      }
     for (k=1; k<=cptcovprod;k++){ /* For product without age */      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]); */        /* 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][1]]==0){
         if(Dummy[Tvard[k][2]==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,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
         }else{          }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,k)] * Tqresult[nres][k];
         }          }
       }else{        }else{
         if(Dummy[Tvard[k][2]==0]){          if(Dummy[Tvard[k][2]]==0){
           cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
         }else{          }else{
           cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];            cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
Line 3121  double **pmij(double **ps, double *cov, Line 3156  double **pmij(double **ps, double *cov,
     ps[i][i]=1./(s1+1.);      ps[i][i]=1./(s1+1.);
     /* Computing other pijs */      /* Computing other pijs */
     for(j=1; j<i; j++)      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++)      for(j=i+1; j<=nlstate+ndeath; j++)
       ps[i][j]= exp(ps[i][j])*ps[i][i];        ps[i][j]= exp(ps[i][j])*ps[i][i];
     /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */      /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
Line 3179  double **pmij(double **ps, double *cov, Line 3214  double **pmij(double **ps, double *cov,
   /* dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); */    /* dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
   
   /* P_x */    /* 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 */    /* outputs pmmij which is a stochastic matrix in row */
   
   /* Diag(w_x) */    /* Diag(w_x) */
Line 3425  double ***hpxij(double ***po, int nhstep Line 3460  double ***hpxij(double ***po, int nhstep
       for (k=1; k<=cptcovprod;k++){ /*  For product without age */        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]); */          /* 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)]; */          /* 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][1]]==0){
           if(Dummy[Tvard[k][2]==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,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
           }else{            }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,k)] * Tqresult[nres][k];
           }            }
         }else{          }else{
           if(Dummy[Tvard[k][2]==0]){            if(Dummy[Tvard[k][2]]==0){
             cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];              cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
           }else{            }else{
             cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];              cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
Line 3524  double ***hbxij(double ***po, int nhstep Line 3559  double ***hbxij(double ***po, int nhstep
       cov[2]=agexact;        cov[2]=agexact;
       if(nagesqr==1)        if(nagesqr==1)
         cov[3]= agexact*agexact;          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,k)]; */
       /* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[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)); */          /* 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 */        for (k=1; k<=nsq;k++) { /* For single varying covariates only */
Line 3546  double ***hbxij(double ***po, int nhstep Line 3581  double ***hbxij(double ***po, int nhstep
       }        }
       for (k=1; k<=cptcovprod;k++){ /* Useless because included in cptcovn */        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,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("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
       /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/        /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
Line 3555  double ***hbxij(double ***po, int nhstep Line 3603  double ***hbxij(double ***po, int nhstep
       /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */        /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */
       /*                                                 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */        /*                                                 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */
       out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\        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){ */        /* if((int)age == 70){ */
       /*        printf(" Backward hbxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */        /*        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++) { */        /*        for(i=1; i<=nlstate+ndeath; i++) { */
Line 3645  double func( double *x) Line 3693  double func( double *x)
         /* # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi */          /* # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi */
         /*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */          /*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
         /*  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 */          /*  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)  */          /* TvarFind;  TvarFind[1]=6,  TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod)  */
         cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (TvarFind[1]=6)*/          cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (TvarFind[1]=6)*/
         /* V1*V2 (7)  TvarFind[2]=7, TvarFind[3]=9 */          /* V1*V2 (7)  TvarFind[2]=7, TvarFind[3]=9 */
       }        }
Line 5861  void  concatwav(int wav[], int **dh, int Line 5909  void  concatwav(int wav[], int **dh, int
   
 {  {
   /* Health expectancies, no variances */    /* 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 i, j, nhstepm, hstepm, h, nstepm;
   int nhstepma, nstepma; /* Decreasing with age */    int nhstepma, nstepma; /* Decreasing with age */
   double age, agelim, hf;    double age, agelim, hf;
Line 6118  void  concatwav(int wav[], int **dh, int Line 6168  void  concatwav(int wav[], int **dh, int
             varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;              varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;
       }        }
     }      }
     if((int)age ==50){      /* if((int)age ==50){ */
       printf(" age=%d cij=%d nres=%d varhe[%d][%d]=%f ",(int)age, cij, nres, 1,2,varhe[1][2]);      /*   printf(" age=%d cij=%d nres=%d varhe[%d][%d]=%f ",(int)age, cij, nres, 1,2,varhe[1][2]); */
     }      /* } */
     /* Computing expectancies */      /* Computing expectancies */
     hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij,nres);        hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij,nres);  
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
Line 6783  void varprob(char optionfilefiname[], do Line 6833  void varprob(char optionfilefiname[], do
    int k2, l2, j1,  z1;     int k2, l2, j1,  z1;
    int k=0, l;     int k=0, l;
    int first=1, first1, first2;     int first=1, first1, first2;
      int nres=0; /* New */
    double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp;     double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp;
    double **dnewm,**doldm;     double **dnewm,**doldm;
    double *xp;     double *xp;
Line 6871  To be simple, these graphs help to under Line 6922  To be simple, these graphs help to under
    if (cptcovn<1) {tj=1;ncodemax[1]=1;}     if (cptcovn<1) {tj=1;ncodemax[1]=1;}
    j1=0;     j1=0;
    for(j1=1; j1<=tj;j1++){  /* For each valid combination of covariates or only once*/     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) {       if  (cptcovn>0) {
        fprintf(ficresprob, "\n#********** Variable ");          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,z1)]);
Line 6906  To be simple, these graphs help to under Line 6959  To be simple, these graphs help to under
        cov[2]=age;         cov[2]=age;
        if(nagesqr==1)         if(nagesqr==1)
          cov[3]= age*age;           cov[3]= age*age;
        for (k=1; k<=cptcovn;k++) {         /* for (k=1; k<=cptcovn;k++) { */
          cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,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,k)];
          /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4           /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4
                                                                     * 1  1 1 1 1                                                                      * 1  1 1 1 1
                                                                     * 2  2 1 1 1                                                                      * 2  2 1 1 1
Line 6918  To be simple, these graphs help to under Line 6974  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 */         /* 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] */         /* ) 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]]=cov[2+Tage[k]]*cov[2]; */
        for (k=1; k<=cptcovage;k++)         for (k=1; k<=cptcovage;k++){  /* For product with age */
          cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];           if(Dummy[Tage[k]]==2){ /* dummy with age */
        for (k=1; k<=cptcovprod;k++)             cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(j1,k)]*cov[2];
          cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];             /* 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,k)] * nbcode[Tvard[k][2]][codtabm(j1,k)];
                /* 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,k)] * 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,k)] * 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(theta=1; theta <=npar; theta++){
          for(i=1; i<=npar; i++)           for(i=1; i<=npar; i++)
            xp[i] = x[i] + (i==theta ?delti[theta]:(double)0);             xp[i] = x[i] + (i==theta ?delti[theta]:(double)0);
Line 7108  To be simple, these graphs help to under Line 7191  To be simple, these graphs help to under
          } /* k12 */           } /* k12 */
        } /*l1 */         } /*l1 */
      }/* k1 */       }/* k1 */
      } /* loop on nres */
    }  /* loop on combination of covariates j1 */     }  /* loop on combination of covariates j1 */
    free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);     free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
    free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);     free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
Line 7237  void printinghtml(char fileresu[], char Line 7321  void printinghtml(char fileresu[], char
       }        }
                 
        /* if(nqfveff+nqtveff 0) */ /* Test to be done */         /* if(nqfveff+nqtveff 0) */ /* Test to be done */
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtm," (model=%s) ************\n<hr size=\"2\" color=\"#EC5E5E\">",model);
        if(invalidvarcomb[k1]){         if(invalidvarcomb[k1]){
          fprintf(fichtm,"\n<h3>Combination (%d) ignored because no cases </h3>\n",k1);            fprintf(fichtm,"\n<h3>Combination (%d) ignored because no cases </h3>\n",k1); 
          printf("\nCombination (%d) ignored because no cases \n",k1);            printf("\nCombination (%d) ignored because no cases \n",k1); 
Line 7258  divided by h: <sub>h</sub>P<sub>ij</sub> Line 7342  divided by h: <sub>h</sub>P<sub>ij</sub>
 <img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);   <img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); 
      /* Survival functions (period) in state j */       /* Survival functions (period) in state j */
      for(cpt=1; cpt<=nlstate;cpt++){       for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \         fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
 <img src=\"%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," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
          fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
      }       }
      /* State specific survival functions (period) */       /* State specific survival functions (period) */
      for(cpt=1; cpt<=nlstate;cpt++){       for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Survival functions in state %d and in any other live state (total).\         fprintf(fichtm,"<br>\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.       \   And probability to be observed in various states (up to %d) being in state %d at different ages.       \
  <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> <img src=\"%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);   <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
          fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
          fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
      }       }
      /* Period (forward stable) prevalence in each health state */       /* Period (forward stable) prevalence in each health state */
      for(cpt=1; cpt<=nlstate;cpt++){       for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\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. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \         fprintf(fichtm,"<br>\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. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
 <img src=\"%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," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"P_"),subdirf2(optionfilefiname,"P_"));
         fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
      }       }
      if(prevbcast==1){       if(prevbcast==1){
        /* Backward prevalence in each health state */         /* Backward prevalence in each health state */
Line 7424  See page 'Matrix of variance-covariance Line 7512  See page 'Matrix of variance-covariance
         fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);          fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
       }        }
   
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtm," (model=%s) ************\n<hr size=\"2\" color=\"#EC5E5E\">",model);
   
        if(invalidvarcomb[k1]){         if(invalidvarcomb[k1]){
          fprintf(fichtm,"\n<h4>Combination (%d) ignored because no cases </h4>\n",k1);            fprintf(fichtm,"\n<h4>Combination (%d) ignored because no cases </h4>\n",k1); 
Line 7560  void printinggnuplot(char fileresu[], ch Line 7648  void printinggnuplot(char fileresu[], ch
         fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);          fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
         fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);          fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);
         /* fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); */          /* fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); */
         fprintf(ficgp,"set title \"Alive state %d %s\" font \"Helvetica,12\"\n",cpt,gplotlabel);          fprintf(ficgp,"set title \"Alive state %d %s model=%s\" font \"Helvetica,12\"\n",cpt,gplotlabel,model);
         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres);          fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres);
         /* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */          /* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */
       /* k1-1 error should be nres-1*/        /* k1-1 error should be nres-1*/
Line 7981  set ter svg size 640, 480\nunset log y\n Line 8069  set ter svg size 640, 480\nunset log y\n
             fprintf(ficgp,", '' ");              fprintf(ficgp,", '' ");
           /* l=(nlstate+ndeath)*(i-1)+1; */            /* l=(nlstate+ndeath)*(i-1)+1; */
           l=(nlstate+ndeath)*(cpt-1)+1; /* fixed for i; cpt=1 1, cpt=2 1+ nlstate+ndeath, 1+2*(nlstate+ndeath) */            l=(nlstate+ndeath)*(cpt-1)+1; /* fixed for i; cpt=1 1, cpt=2 1+ nlstate+ndeath, 1+2*(nlstate+ndeath) */
           /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */            /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */
           /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */            /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */
           fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+i-1); /* To be verified */            fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+i-1); /* To be verified */
           /* for (j=2; j<= nlstate ; j ++) */            /* for (j=2; j<= nlstate ; j ++) */
           /*    fprintf(ficgp,"+$%d",k+l+j-1); */            /*    fprintf(ficgp,"+$%d",k+l+j-1); */
Line 8331  set ter svg size 640, 480\nunset log y\n Line 8419  set ter svg size 640, 480\nunset log y\n
             /* for(j=3; j <=ncovmodel-nagesqr; j++) { */              /* for(j=3; j <=ncovmodel-nagesqr; j++) { */
             for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */              for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */
               /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */                /* 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 */                switch(Typevar[j]){
                 if(j==Tage[ij]) { /* Product by age  To be looked at!!*/                case 1:
                   if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */                  if(cptcovage >0){ /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
                     if(DummyV[j]==0){                    if(j==Tage[ij]) { /* Product by age  To be looked at!!*//* Bug valgrind */
                       fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);;                      if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
                     }else{ /* quantitative */                        if(DummyV[j]==0){/* Bug valgrind */
                       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+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);;
                       /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,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*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
                         }
                         ij++;
                     }                      }
                     ij++;  
                   }                    }
                 }                   }
               }else if(cptcovprod >0){                  break;
                 if(j==Tprod[ijp]) { /* */                 case 2:
                   /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */                  if(cptcovprod >0){
                   if(ijp <=cptcovprod) { /* Product */                    if(j==Tprod[ijp]) { /* */ 
                     if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */                      /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
                       if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */                      if(ijp <=cptcovprod) { /* Product */
                         /* 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)]); */                        if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */
                         fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]);                          if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */
                       }else{ /* Vn is dummy and Vm is quanti */                            /* 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*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[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]]);
                         fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[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]]); */
                     }else{ /* Vn*Vm Vn is quanti */                            fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
                       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{ /* Vn*Vm Vn is quanti */
                       }else{ /* Both quanti */                          if(DummyV[Tvard[ijp][2]]==0){
                         fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);                            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 */
                   }                  }
                 } /* end Tprod */                  break;
               } else{  /* simple covariate */                case 0:
                   /* simple covariate */
                 /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(k1,j)]); /\* Valgrind bug nbcode *\/ */                  /* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(k1,j)]); /\* Valgrind bug nbcode *\/ */
                 if(Dummy[j]==0){                  if(Dummy[j]==0){
                   fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]); /*  */                    fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]); /*  */
Line 8373  set ter svg size 640, 480\nunset log y\n Line 8469  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*%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])]); */                    /* 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 */              } /* end j */
           }else{            }else{ /* k=k2 */
             i=i-ncovmodel;              if(ng !=1 ){ /* For logit formula of log p11 is more difficult to get */
             if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */                fprintf(ficgp," (1.");i=i-ncovmodel;
               fprintf(ficgp," (1.");              }else
                 i=i-ncovmodel;
           }            }
                       
           if(ng != 1){            if(ng != 1){
Line 8391  set ter svg size 640, 480\nunset log y\n Line 8492  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);                  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;                ij=1;
               for(j=3; j <=ncovmodel-nagesqr; j++){                ijp=1;
                  if(cptcovage >0){                 /* for(j=3; j <=ncovmodel-nagesqr; j++){ */
                    if((j-2)==Tage[ij]) { /* Bug valgrind */                for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */
                      if(ij <=cptcovage) { /* Bug valgrind */                  switch(Typevar[j]){
                        fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);                  case 1:
                        /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */                    if(cptcovage >0){ 
                        ij++;                      if(j==Tage[ij]) { /* Bug valgrind */
                      }                        if(ij <=cptcovage) { /* Bug valgrind */
                    }                          if(DummyV[j]==0){/* Bug valgrind */
                  }else                            /* 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",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);/* Valgrind bug nbcode */                            /* 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,")");                fprintf(ficgp,")");
             }              }
Line 8410  set ter svg size 640, 480\nunset log y\n Line 8572  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);                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 */              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);                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 */              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);                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);
           }            }
Line 9524  int readdata(char datafile[], int firsto Line 9686  int readdata(char datafile[], int firsto
   
   DummyV=ivector(1,NCOVMAX); /* 1 to 3 */    DummyV=ivector(1,NCOVMAX); /* 1 to 3 */
   FixedV=ivector(1,NCOVMAX); /* 1 to 3 */    FixedV=ivector(1,NCOVMAX); /* 1 to 3 */
     for(v=1;v<NCOVMAX;v++){
       DummyV[v]=0;
       FixedV[v]=0;
     }
   
   for(v=1; v <=ncovcol;v++){    for(v=1; v <=ncovcol;v++){
     DummyV[v]=0;      DummyV[v]=0;
Line 11163  int hPijx(double *p, int bage, int fage) Line 11329  int hPijx(double *p, int bage, int fage)
   
         /* oldm=oldms;savm=savms; */          /* oldm=oldms;savm=savms; */
         /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);   */          /* 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); */          /* 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=");          fprintf(ficrespijb,"# Cov Agex agex-h hbijx with i,j=");
         for(i=1; i<=nlstate;i++)          for(i=1; i<=nlstate;i++)
Line 11176  int hPijx(double *p, int bage, int fage) Line 11342  int hPijx(double *p, int bage, int fage)
           /* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */            /* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */
           for(i=1; i<=nlstate;i++)            for(i=1; i<=nlstate;i++)
             for(j=1; j<=nlstate+ndeath;j++)              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");            fprintf(ficrespijb,"\n");
         }          }
         free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);          free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
Line 11861  Please run with mle=-1 to get a correct Line 12027  Please run with mle=-1 to get a correct
   }    }
   mint=matrix(1,maxwav,firstobs,lastobs);    mint=matrix(1,maxwav,firstobs,lastobs);
   anint=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);    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 */    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 */    ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */
Line 11925  Please run with mle=-1 to get a correct Line 12092  Please run with mle=-1 to get a correct
   Tage=ivector(1,NCOVMAX); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age    Tage=ivector(1,NCOVMAX); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age
                          4 covariates (3 plus signs)                           4 covariates (3 plus signs)
                          Tage[1=V3*age]= 4; Tage[2=age*V4] = 3                           Tage[1=V3*age]= 4; Tage[2=age*V4] = 3
                       */                               */  
     for(i=1;i<NCOVMAX;i++)
       Tage[i]=0;
   Tmodelind=ivector(1,NCOVMAX);/** gives the k model position of an    Tmodelind=ivector(1,NCOVMAX);/** gives the k model position of an
                                 * individual dummy, fixed or varying:                                  * individual dummy, fixed or varying:
                                 * Tmodelind[Tvaraff[3]]=9,Tvaraff[1]@9={4,                                  * Tmodelind[Tvaraff[3]]=9,Tvaraff[1]@9={4,
Line 12029  Please run with mle=-1 to get a correct Line 12198  Please run with mle=-1 to get a correct
            * For k=4 covariates, h goes from 1 to m=2**k             * For k=4 covariates, h goes from 1 to m=2**k
            * codtabm(h,k)=  (1 & (h-1) >> (k-1)) + 1;             * codtabm(h,k)=  (1 & (h-1) >> (k-1)) + 1;
            * #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1             * #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
            *     h\k   1     2     3     4             *     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             *     1 i=1 1 i=1 1 i=1 1 i=1 1   *     0     0  0  0  0 
            *     2     2     1     1     1             *     2     2     1     1     1   *     1     0  0  0  1 
            *     3 i=2 1     2     1     1             *     3 i=2 1     2     1     1   *     2     0  0  1  0 
            *     4     2     2     1     1             *     4     2     2     1     1   *     3     0  0  1  1 
            *     5 i=3 1 i=2 1     2     1             *     5 i=3 1 i=2 1     2     1   *     4     0  1  0  0 
            *     6     2     1     2     1             *     6     2     1     2     1   *     5     0  1  0  1 
            *     7 i=4 1     2     2     1             *     7 i=4 1     2     2     1   *     6     0  1  1  0 
            *     8     2     2     2     1             *     8     2     2     2     1   *     7     0  1  1  1 
            *     9 i=5 1 i=3 1 i=2 1     2             *     9 i=5 1 i=3 1 i=2 1     2   *     8     1  0  0  0 
            *    10     2     1     1     2             *    10     2     1     1     2   *     9     1  0  0  1 
            *    11 i=6 1     2     1     2             *    11 i=6 1     2     1     2   *    10     1  0  1  0 
            *    12     2     2     1     2             *    12     2     2     1     2   *    11     1  0  1  1 
            *    13 i=7 1 i=4 1     2     2                 *    13 i=7 1 i=4 1     2     2   *    12     1  1  0  0  
            *    14     2     1     2     2             *    14     2     1     2     2   *    13     1  1  0  1 
            *    15 i=8 1     2     2     2             *    15 i=8 1     2     2     2   *    14     1  1  1  0 
            *    16     2     2     2     2             *    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? */    /* 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       /* from h=5 and m, we get then number of covariates k=log(m)/log(2)=4
      * and the value of each covariate?       * and the value of each covariate?
Line 12138  Title=%s <br>Datafile=%s Firstpass=%d La Line 12307  Title=%s <br>Datafile=%s Firstpass=%d La
           optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);            optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }    }
   
   fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C)  2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-2013-2016-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br>  \    fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C)  2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-2013-2016-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br>  \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
 <font size=\"2\">IMaCh-%s <br> %s</font> \  <font size=\"2\">IMaCh-%s <br> %s</font> \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
Line 12548  Please run with mle=-1 to get a correct Line 12717  Please run with mle=-1 to get a correct
       hesscov(matcov, hess, p, npar, delti, ftolhess, func);        hesscov(matcov, hess, p, npar, delti, ftolhess, func);
       printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");        printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
       fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n  It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");        fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n  It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
         fprintf(fichtm, "\n<p>The Wald test results are output only if the maximimzation of the Likelihood is performed (mle=1)\n</br>Parameters, Wald tests and Wald-based confidence intervals\n</br> W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n</br> And Wald-based confidence intervals plus and minus 1.96 * W \n </br> It might be better to visualize the covariance matrix. See the page '<a href=\"%s\">Matrix of variance-covariance of one-step probabilities and its graphs</a>'.\n</br>",optionfilehtmcov);
       fprintf(fichtm,"\n<table style=\"text-align:center; border: 1px solid\">");        fprintf(fichtm,"\n<table style=\"text-align:center; border: 1px solid\">");
       fprintf(fichtm, "\n<tr><th>Model=</th><th>1</th><th>+ age</th>");        fprintf(fichtm, "\n<tr><th>Model=</th><th>1</th><th>+ age</th>");
       if(nagesqr==1){        if(nagesqr==1){
Line 12581  Please run with mle=-1 to get a correct Line 12751  Please run with mle=-1 to get a correct
               printf("%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));                printf("%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
               fprintf(ficlog,"%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));                fprintf(ficlog,"%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
               if(fabs(wald) > 1.96){                if(fabs(wald) > 1.96){
                 fprintf(fichtm, "<td><b>%12.7f</b> (%12.7f)</br>",p[jk],sqrt(matcov[jk][jk]));                  fprintf(fichtm, "<td><b>%12.7f</b></br> (%12.7f)</br>",p[jk],sqrt(matcov[jk][jk]));
                 fprintf(fichtm,"<b>W=%8.3f</b></br>",wald);  
               }else{                }else{
                 fprintf(fichtm, "<td>%12.7f (%12.7f)</br>",p[jk],sqrt(matcov[jk][jk]));                  fprintf(fichtm, "<td>%12.7f (%12.7f)</br>",p[jk],sqrt(matcov[jk][jk]));
                 fprintf(fichtm,"W=%8.3f</br>",wald);  
               }                }
                 fprintf(fichtm,"W=%8.3f</br>",wald);
               fprintf(fichtm,"[%12.7f;%12.7f]</br></td>", p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));                fprintf(fichtm,"[%12.7f;%12.7f]</br></td>", p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
               jk++;                 jk++; 
             }              }
Line 13188  Please run with mle=-1 to get a correct Line 13357  Please run with mle=-1 to get a correct
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */      for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
       if(i1 != 1 && TKresult[nres]!= k)        if(i1 != 1 && TKresult[nres]!= k)
         continue;          continue;
       printf("\n#****** Result for:");        printf("\n# model %s \n#****** Result for:", model);
       fprintf(ficrest,"\n#****** Result for:");        fprintf(ficrest,"\n# model %s \n#****** Result for:", model);
       fprintf(ficlog,"\n#****** Result for:");        fprintf(ficlog,"\n# model %s \n#****** Result for:", model);
       for(j=1;j<=cptcoveff;j++){         for(j=1;j<=cptcoveff;j++){ 
         printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,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(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
Removed from v.1.319  
changed lines
  Added in v.1.329

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>