imach/src/imach.c - diff

Return to imach.c CVS log

Up to [Local Repository] / imach / src

Diff for /imach/src/imach.c between versions 1.339 and 1.352

-version 1.339, 2022/09/09 17:55:22
+version 1.352, 2023/04/29 10:46:21
  Line 1
  /* $Id$
    $State$
    $Log$
+   Revision 1.352  2023/04/29 10:46:21  brouard
+   *** empty log message ***
+   Revision 1.351  2023/04/29 10:43:47  brouard
+   Summary: 099r45
+   Revision 1.350  2023/04/24 11:38:06  brouard
+   *** empty log message ***
+   Revision 1.349  2023/01/31 09:19:37  brouard
+   Summary: Improvements in models with age*Vn*Vm
+   Revision 1.347  2022/09/18 14:36:44  brouard
+   Summary: version 0.99r42
+   Revision 1.346  2022/09/16 13:52:36  brouard
+   * src/imach.c (Module): 0.99r41 Was an error when product of timevarying and fixed. Using FixedV[of name] now. Thank you  Feinuo
+   Revision 1.345  2022/09/16 13:40:11  brouard
+   Summary: Version 0.99r41
+   * imach.c (Module): 0.99r41 Was an error when product of timevarying and fixed. Using FixedV[of name] now. Thank you  Feinuo
+   Revision 1.344  2022/09/14 19:33:30  brouard
+   Summary: version 0.99r40
+   * imach.c (Module): Fixing names of variables in T_ (thanks to Feinuo)
+   Revision 1.343  2022/09/14 14:22:16  brouard
+   Summary: version 0.99r39
+   * imach.c (Module): Version 0.99r39 with colored dummy covariates
+   (fixed or time varying), using new last columns of
+   ILK_parameter.txt file.
+   Revision 1.342  2022/09/11 19:54:09  brouard
+   Summary: 0.99r38
+   * imach.c (Module): Adding timevarying products of any kinds,
+   should work before shifting cotvar from ncovcol+nqv columns in
+   order to have a correspondance between the column of cotvar and
+   the id of column.
+   (Module): Some cleaning and adding covariates in ILK.txt
+   Revision 1.341  2022/09/11 07:58:42  brouard
+   Summary: Version 0.99r38
+   After adding change in cotvar.
+   Revision 1.340  2022/09/11 07:53:11  brouard
+   Summary: Version imach 0.99r37
+   * imach.c (Module): Adding timevarying products of any kinds,
+   should work before shifting cotvar from ncovcol+nqv columns in
+   order to have a correspondance between the column of cotvar and
+   the id of column.
    Revision 1.339  2022/09/09 17:55:22  brouard
    Summary: version 0.99r37
- Line 1269  typedef struct {
+ Line 1326  typedef struct {
  /* #include <libintl.h> */
  /* #define _(String) gettext (String) */
- #define MAXLINE 2048 /* Was 256 and 1024. Overflow with 312 with 2 states and 4 covariates. Should be ok */
+ #define MAXLINE 16384 /* Was 256 and 1024 and 2048. Overflow with 312 with 2 states and 4 covariates. Should be ok */
  #define GNUPLOTPROGRAM "gnuplot"
+ #define GNUPLOTVERSION 5.1
+ double gnuplotversion=GNUPLOTVERSION;
  /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
  #define FILENAMELENGTH 256
  #define GLOCK_ERROR_NOPATH              -1      /* empty path */
  #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */
- #define MAXPARM 128 /**< Maximum number of parameters for the optimization */
+ #define MAXPARM 216 /**< Maximum number of parameters for the optimization was 128 */
  #define NPARMAX 64 /**< (nlstate+ndeath-1)*nlstate*ncovmodel */
  #define NINTERVMAX 8
- Line 1312  typedef struct {
+ Line 1371  typedef struct {
  /* $State$ */
  #include "version.h"
  char version[]=__IMACH_VERSION__;
- char copyright[]="September 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[]="April 2023,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 strstart[80];
  char optionfilext[10], optionfilefiname[FILENAMELENGTH];
  int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
+ int debugILK=0; /* debugILK is set by a #d in a comment line */
  int nagesqr=0, nforce=0; /* nagesqr=1 if model is including age*age, number of forces */
  /* 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 */
- Line 1325  int cptcovt=0; /**< cptcovt: total numbe
+ Line 1385  int cptcovt=0; /**< cptcovt: total numbe
  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 cptcovprodage=0; /**< Number of fixed covariates with age: V3*age or V2*V3*age 1 */
+ int cptcovprodvage=0; /**< Number of varying covariates with age: V7*age or V7*V6*age */
+ int cptcovdageprod=0; /**< Number of doubleproducts with age, since 0.99r44 only: age*Vn*Vm for gnuplot printing*/
  int cptcovprodnoage=0; /**< Number of covariate products without age */
  int cptcoveff=0; /* Total number of single dummy covariates (fixed or time varying) 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 ncovvt=0; /* Total number of effective (wave) varying covariates (dummy or quantitative or products [without age]) in the model */
- int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (dummy of quantitative) in the model */
+ int ncovvta=0; /*  +age*V6 + age*V7+ age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 Total number of expandend products [with age]) in the model */
+ int ncovta=0; /*age*V3*V2 +age*V2+agev3+ageV4  +age*V6 + age*V7+ age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 Total number of expandend products [with age]) in the model */
+ int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (single or product, dummy or quantitative) in the model */
+ int ncovva=0; /* +age*V6 + age*V7+ge*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 Total number of effective (wave and stepm) varying with age covariates (single or product, dummy or quantitative) in the model */
  int nsd=0; /**< Total number of single dummy variables (output) */
  int nsq=0; /**< Total number of single quantitative variables (output) */
  int ncoveff=0; /* Total number of effective fixed dummy covariates in the model */
- Line 1417  extern time_t time();
+ Line 1483  extern time_t time();
  struct tm start_time, end_time, curr_time, last_time, forecast_time;
  time_t  rstart_time, rend_time, rcurr_time, rlast_time, rforecast_time; /* raw time */
+ time_t   rlast_btime; /* raw time */
  struct tm tm;
  char strcurr[80], strfor[80];
- Line 1495  double  **covar; /**< covar[j,i], value
+ Line 1562  double  **covar; /**< covar[j,i], value
                    * covar=matrix(0,NCOVMAX,1,n);
                    * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */
  double **coqvar; /* Fixed quantitative covariate nqv */
- double ***cotvar; /* Time varying covariate ntv */
+ double ***cotvar; /* Time varying covariate start at ncovcol + nqv + (1 to ntv) */
  double ***cotqvar; /* Time varying quantitative covariate itqv */
  double  idx;
  int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
- Line 1507  int **nbcode, *Tvar; /**< model=V2 => Tv
+ Line 1574  int **nbcode, *Tvar; /**< model=V2 => Tv
         *  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
+        *       cotvar[wav][ncovcol+nqv+ iv(1 to nqtv)][i]= [(1 to nqtv)][i]=(V12) quanti
         *       cotvar[wav][iv(1 to ntv)][i]= [1][1][i]=(V9) dummies at wav 1
         *       cotvar[wav][iv(1 to ntv)][i]= [1][2][i]=(V10) dummies at wav 1
         *       covar[Vk,i], value of the Vkth fixed covariate dummy or quanti for individual i:
- Line 1520  int **nbcode, *Tvar; /**< model=V2 => Tv
+ Line 1587  int **nbcode, *Tvar; /**< model=V2 => Tv
    # States 1=Coresidence, 2 Living alone, 3 Institution
    # 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+V4*V3*age */
- /*    k        1  2   3   4     5    6    7     8    9 */
+ /*    kmodel  1  2   3    4     5     6    7     8     9    10 */
- /*Typevar[k]=  0  0   0   2     1    0    2     1    0 *//*0 for simple covariate (dummy, quantitative,*/
+ /*Typevar[k]=  0  0   0   2     1    0    2     1     0    3 *//*0 for simple covariate (dummy, quantitative,*/
-                                                          /* fixed or varying), 1 for age product, 2 for*/
+                                                                /* fixed or varying), 1 for age product, 2 for*/
-                                                          /* product */
+                                                                /* product without age, 3 for age and double product   */
- /*Dummy[k]=    1  0   0   1     3    1    1     2    0 *//*Dummy[k] 0=dummy (0 1), 1 quantitative */
+ /*Dummy[k]=    1  0   0   1     3    1    1     2     0     3  *//*Dummy[k] 0=dummy (0 1), 1 quantitative */
-                                                          /*(single or product without age), 2 dummy*/
+                                                                 /*(single or product without age), 2 dummy*/
-                                                          /* with age product, 3 quant with age product*/
+                                                                /* with age product, 3 quant with age product*/
- /*Tvar[k]=     5  4   3   6     5    2    7     1    1 */
+ /*Tvar[k]=     5  4   3   6     5    2    7     1     1     6 */
- /*    nsd         1   2                              3 */ /* Counting single dummies covar fixed or tv */
+ /*    nsd         1   2                               3 */ /* Counting single dummies covar fixed or tv */
- /*TnsdVar[Tvar]   1   2                              3 */
+ /*TnsdVar[Tvar]   1   2                               3 */
- /*Tvaraff[nsd]     4   3                              1 */ /* ID of single dummy cova fixed or timevary*/
+ /*Tvaraff[nsd]    4   3                               1 */ /* ID of single dummy cova fixed or timevary*/
- /*TvarsD[nsd]     4   3                              1 */ /* ID of single dummy cova fixed or timevary*/
+ /*TvarsD[nsd]     4   3                               1 */ /* ID of single dummy cova fixed or timevary*/
- /*TvarsDind[nsd]  2   3                              9 */ /* position K of single dummy cova */
+ /*TvarsDind[nsd]  2   3                               9 */ /* position K of single dummy cova */
- /*    nsq      1                     2                 */ /* Counting single quantit tv */
+ /*    nsq      1                     2                  */ /* Counting single quantit tv */
- /* TvarsQ[k]   5                     2                 */ /* Number of single quantitative cova */
+ /* TvarsQ[k]   5                     2                  */ /* Number of single quantitative cova */
- /* TvarsQind   1                     6                 */ /* position K of single quantitative cova */
+ /* TvarsQind   1                     6                  */ /* position K of single quantitative cova */
- /* Tprod[i]=k             1               2            */ /* Position in model of the ith prod without age */
+ /* Tprod[i]=k             1               2             */ /* Position in model of the ith prod without age */
- /* cptcovage                    1               2      */ /* Counting cov*age in the model equation */
+ /* cptcovage                    1               2         3 */ /* Counting cov*age in the model equation */
- /* Tage[cptcovage]=k            5               8      */ /* Position in the model of ith cov*age */
+ /* Tage[cptcovage]=k            5               8         10 */ /* 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 */
+ /* model="V2+V3+V4+V6+V7+V6*V2+V7*V2+V6*V3+V7*V3+V6*V4+V7*V4+age*V2+age*V3+age*V4+age*V6+age*V7+age*V6*V2+age*V6*V3+age*V7*V3+age*V6*V4+age*V7*V4\r"*/
+ /*  p Tvard[1][1]@21 = {6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0}*/
+ /*  p Tvard[2][1]@21 = {7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0 <repeats 11 times>}
+ /* p Tvardk[1][1]@24 = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0}*/
+ /* p Tvardk[1][1]@22 = {0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0} */
+ /* 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)  */
- Line 1590  int *TvarVQ; /* TvarVQ[1]=V5 in V5+V4+V3
+ Line 1662  int *TvarVQ; /* TvarVQ[1]=V5 in V5+V4+V3
  int *TvarVQind; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */
  int *TvarVV; /* We count ncovvt time varying covariates (single or products without age) and put their name into TvarVV */
  int *TvarVVind; /* We count ncovvt time varying covariates (single or products without age) and put their name into TvarVV */
+ int *TvarVVA; /* We count ncovvt time varying covariates (single or products with age) and put their name into TvarVVA */
+ int *TvarVVAind; /* We count ncovvt time varying covariates (single or products without age) and put their name into TvarVV */
+ int *TvarAVVA; /* We count ALL ncovta time varying covariates (single or products with age) and put their name into TvarVVA */
+ int *TvarAVVAind; /* We count ALL ncovta time varying covariates (single or products without age) and put their name into TvarVV */
        /*#  ID           V1     V2          weight               birth   death   1st    s1      V3      V4      V5       2nd  s2 */
-       /* model V1+V3+age*V1+age*V3+V1*V3 */
+       /* model V1+V3+age*V1+age*V3+V1*V3 + V1*V3*age */
-       /*  Tvar={1, 3, 1, 3, 6}, the 6 comes from the fact that there are already V1, V2, V3, V4, V5 native covariates */
+       /*  Tvar={1, 3, 1, 3, 6, 6}, the 6 comes from the fact that there are already V1, V2, V3, V4, V5 native covariates */
-       /* TvarVV={3,1,3}, for V3 and then the product V1*V3 is decomposed into V1 and V3 */
+       /* TvarVV={3,1,3,1,3}, for V3 and then the product V1*V3 is decomposed into V1 and V3 */
-       /* TvarVVind={2,5,5}, for V3 and then the product V1*V3 is decomposed into V1 and V3 */
+       /* TvarVVind={2,5,5,6,6}, for V3 and then the product V1*V3 is decomposed into V1 and V3 and V1*V3*age into 6,6 */
  int *Tvarsel; /**< Selected covariates for output */
  double *Tvalsel; /**< Selected modality value of covariate for output */
- int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product */
+ int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product, 3 age*Vn*Vm */
  int *Fixed; /** Fixed[k] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */
  int *Dummy; /** Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */
  int *DummyV; /** Dummy[v] 0=dummy (0 1), 1 quantitative */
- Line 1734  char *trimbb(char *out, char *in)
+ Line 1810  char *trimbb(char *out, char *in)
    return s;
  }
+ char *trimbtab(char *out, char *in)
+ { /* Trim  blanks or tabs in line but keeps first blanks if line starts with blanks */
+   char *s;
+   s=out;
+   while (*in != '\0'){
+     while( (*in == ' ' || *in == '\t')){ /* && *(in+1) != '\0'){*/
+       in++;
+     }
+     *out++ = *in++;
+   }
+   *out='\0';
+   return s;
+ }
  /* char *substrchaine(char *out, char *in, char *chain) */
  /* { */
  /*   /\* Substract chain 'chain' from 'in', return and output 'out' *\/ */
- Line 1760  char *trimbb(char *out, char *in)
+ Line 1850  char *trimbb(char *out, char *in)
  char *substrchaine(char *out, char *in, char *chain)
  {
    /* Substract chain 'chain' from 'in', return and output 'out' */
-   /* in="V1+V1*age+age*age+V2", chain="age*age" */
+   /* in="V1+V1*age+age*age+V2", chain="+age*age" out="V1+V1*age+V2" */
    char *strloc;
-   strcpy (out, in);
+   strcpy (out, in);                   /* out="V1+V1*age+age*age+V2" */
-   strloc = strstr(out, chain); /* strloc points to out at age*age+V2 */
+   strloc = strstr(out, chain); /* strloc points to out at "+age*age+V2"  */
-   printf("Bef strloc=%s chain=%s out=%s \n", strloc, chain, out);
+   printf("Bef strloc=%s chain=%s out=%s \n", strloc, chain, out); /* strloc=+age*age+V2 chain="+age*age", out="V1+V1*age+age*age+V2" */
    if(strloc != NULL){
-     /* will affect out */ /* strloc+strlenc(chain)=+V2 */ /* Will also work in Unicode */
+     /* will affect out */ /* strloc+strlen(chain)=|+V2 = "V1+V1*age+age*age|+V2" */ /* Will also work in Unicodek */
-     memmove(strloc,strloc+strlen(chain), strlen(strloc+strlen(chain))+1);
+     memmove(strloc,strloc+strlen(chain), strlen(strloc+strlen(chain))+1); /* move number of bytes corresponding to the length of "+V2" which is 3, plus one is 4 (including the null)*/
-     /* strcpy (strloc, strloc +strlen(chain));*/
+     /* equivalent to strcpy (strloc, strloc +strlen(chain)) if no overlap; Copies from "+V2" to V1+V1*age+ */
    }
-   printf("Aft strloc=%s chain=%s in=%s out=%s \n", strloc, chain, in, out);
+   printf("Aft strloc=%s chain=%s in=%s out=%s \n", strloc, chain, in, out);  /* strloc=+V2 chain="+age*age", in="V1+V1*age+age*age+V2", out="V1+V1*age+V2" */
    return out;
  }
- Line 1780  char *substrchaine(char *out, char *in,
+ Line 1870  char *substrchaine(char *out, char *in,
  char *cutl(char *blocc, char *alocc, char *in, char occ)
  {
    /* cuts string in into blocc and alocc where blocc ends before FIRST occurence of char 'occ'
-      and alocc starts after first occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2')
+      and alocc starts after first occurence of char 'occ' : ex cutl(blocc,alocc,"abcdef2ghi2j",'2')
       gives alocc="abcdef" and blocc="ghi2j".
       If occ is not found blocc is null and alocc is equal to in. Returns blocc
    */
- Line 1846  int nbocc(char *s, char occ)
+ Line 1936  int nbocc(char *s, char occ)
    return j;
  }
+ int nboccstr(char *textin, char *chain)
+ {
+   /* Counts the number of occurence of "chain"  in string textin */
+   /*  in="+V7*V4+age*V2+age*V3+age*V4"  chain="age" */
+   char *strloc;
+   int i,j=0;
+   i=0;
+   strloc=textin; /* strloc points to "^+V7*V4+age+..." in textin */
+   for(;;) {
+     strloc= strstr(strloc,chain); /* strloc points to first character of chain in textin if found. Example strloc points^ to "+V7*V4+^age" in textin  */
+     if(strloc != NULL){
+       strloc = strloc+strlen(chain); /* strloc points to "+V7*V4+age^" in textin */
+       j++;
+     }else
+       break;
+   }
+   return j;
+ }
  /* void cutv(char *u,char *v, char*t, char occ) */
  /* { */
  /*   /\* cuts string t into u and v where u ends before last occurence of char 'occ'  */
- Line 2525  void powell(double p[], double **xi, int
+ Line 2637  void powell(double p[], double **xi, int
    double fp,fptt;
    double *xits;
    int niterf, itmp;
+   int Bigter=0, nBigterf=1;
    pt=vector(1,n);
    ptt=vector(1,n);
    xit=vector(1,n);
- Line 2538  void powell(double p[], double **xi, int
+ Line 2651  void powell(double p[], double **xi, int
      ibig=0;
      del=0.0;
      rlast_time=rcurr_time;
+     rlast_btime=rcurr_time;
      /* (void) gettimeofday(&curr_time,&tzp); */
      rcurr_time = time(NULL);
      curr_time = *localtime(&rcurr_time);
      /* 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 gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); */
-     printf("\nPowell iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);
+     Bigter=(*iter - *iter % ncovmodel)/ncovmodel +1; /* Big iteration, i.e on ncovmodel cycle */
-     fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
+     printf("\nPowell iter=%d Big Iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,Bigter,*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);
- /*     fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */
+     fprintf(ficlog,"\nPowell iter=%d Big Iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,Bigter,*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
+     fprintf(ficrespow,"%d %d %.12f %d",*iter,Bigter, *fret,curr_time.tm_sec-start_time.tm_sec);
      fp=(*fret); /* From former iteration or initial value */
      for (i=1;i<=n;i++) {
        fprintf(ficrespow," %.12lf", p[i]);
- Line 2567  void powell(double p[], double **xi, int
+ Line 2682  void powell(double p[], double **xi, int
        }else if(Typevar[j]==2) {
          printf("  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
          fprintf(ficlog,"  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+       }else if(Typevar[j]==3) {
+         printf("  +    V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+         fprintf(ficlog,"  +    V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
        }
      }
      printf("\n");
- Line 2597  void powell(double p[], double **xi, int
+ Line 2715  void powell(double p[], double **xi, int
          strcurr[itmp-1]='\0';
        printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
        fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
-       for(niterf=10;niterf<=30;niterf+=10){
+       for(nBigterf=1;nBigterf<=31;nBigterf+=10){
+         niterf=nBigterf*ncovmodel;
+         /* rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time); */
          rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);
          forecast_time = *localtime(&rforecast_time);
          strcpy(strfor,asctime(&forecast_time));
          itmp = strlen(strfor);
          if(strfor[itmp-1]=='\n')
            strfor[itmp-1]='\0';
-         printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
+         printf("   - if your program needs %d BIG iterations (%d iterations) to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",nBigterf, niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
-         fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
+         fprintf(ficlog,"   - if your program needs %d BIG iterations  (%d iterations) to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",nBigterf, niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
        }
      }
      for (i=1;i<=n;i++) { /* For each direction i */
- Line 2906  void powell(double p[], double **xi, int
+ Line 3026  void powell(double p[], double **xi, int
       /* 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(Typevar[k1]==1){ /* A product with age */
+        if(Typevar[k1]==1 || Typevar[k1]==3){ /* A product with age */
           cov[2+nagesqr+k1]=precov[nres][k1]*cov[2];
         }else{
           cov[2+nagesqr+k1]=precov[nres][k1];
- Line 3116  void powell(double p[], double **xi, int
+ Line 3236  void powell(double p[], double **xi, int
        cov[3]= agefin*agefin;;
      }
      for(k1=1;k1<=cptcovt;k1++){ /* loop on model equation (including products) */
-       if(Typevar[k1]==1){ /* A product with age */
+       if(Typevar[k1]==1 || Typevar[k1]==3){ /* A product with age */
          cov[2+nagesqr+k1]=precov[nres][k1]*cov[2];
        }else{
          cov[2+nagesqr+k1]=precov[nres][k1];
- Line 3583  double ***hpxij(double ***po, int nhstep
+ Line 3703  double ***hpxij(double ***po, int nhstep
        /* 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(Typevar[k1]==1){ /* A product with age */
+         if(Typevar[k1]==1 || Typevar[k1]==3){ /* A product with age */
            cov[2+nagesqr+k1]=precov[nres][k1]*cov[2];
          }else{
            cov[2+nagesqr+k1]=precov[nres][k1];
- Line 3769  double ***hbxij(double ***po, int nhstep
+ Line 3889  double ***hbxij(double ***po, int nhstep
        }
        /** New code */
        for(k1=1;k1<=cptcovt;k1++){ /* loop on model equation (including products) */
-         if(Typevar[k1]==1){ /* A product with age */
+         if(Typevar[k1]==1 || Typevar[k1]==3){ /* A product with age */
            cov[2+nagesqr+k1]=precov[nres][k1]*cov[2];
          }else{
            cov[2+nagesqr+k1]=precov[nres][k1];
- Line 3880  double func( double *x)
+ Line 4000  double func( double *x)
    int ioffset=0;
    int ipos=0,iposold=0,ncovv=0;
+   double cotvarv, cotvarvold;
    double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
    double **out;
    double lli; /* Individual log likelihood */
- Line 3911  double func( double *x)
+ Line 4032  double func( double *x)
        */
        ioffset=2+nagesqr ;
     /* Fixed */
-       for (kf=1; kf<=ncovf;kf++){ /* For each fixed covariate dummu or quant or prod */
+       for (kf=1; kf<=ncovf;kf++){ /* For each fixed covariate dummy or quant or prod */
          /* # 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 */
          /*  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 */
- Line 3927  double func( double *x)
+ Line 4048  double func( double *x)
           mw[mi][i] is real wave of the mi th effectve wave */
        /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i];
           s2=s[mw[mi+1][i]][i];
-          And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i]
+          And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i] because now is moved after nvocol+nqv
           But if the variable is not in the model TTvar[iv] is the real variable effective in the model:
           meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]
        */
- Line 3937  double func( double *x)
+ Line 4058  double func( double *x)
          /*   /\* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; but where is the crossproduct? *\/ */
          /*   cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; */
          /* } */
-         for(ncovv=1, ipos=0; ncovv <= ncovvt ; ncovv++){ /* Varying  covariates (single and product but no age )*/
+         for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Varying  covariates (single and product but no age )*/
-           itv=TvarVV[ncovv]; /*  TvarVV={3, 1, 3}  */
+           itv=TvarVV[ncovv]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate */
-           ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] */
+           ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
+           if(FixedV[itv]!=0){ /* Not a fixed covariate */
+             cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i];  /* cotvar[wav][ncovcol+nqv+iv][i] */
+           }else{ /* fixed covariate */
+             cotvarv=covar[itv][i];  /* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model */
+           }
            if(ipos!=iposold){ /* Not a product or first of a product */
-             /* TvarFind={1,0,0,0}  */
+             cotvarvold=cotvarv;
-             if(TvarFind[itv]==0){
+           }else{ /* A second product */
-             cov[ioffset+ipos]= cotvar[mw[mi][i]][ncovv][i];  /* Should be covar if fixed covar[Tvar[TvarFind[itv]]][i]*/
+             cotvarv=cotvarv*cotvarvold;
-             }else{
-               cov[ioffset+ipos]=covar[Tvar[TvarFind[itv]]][i];
-             }
-           }else{
-             if(TvarFind[itv]==0){
-               cov[ioffset+ipos]*= cotvar[mw[mi][i]][ncovv][i];  /* Should be covar if fixed covar[Tvar[TvarFind[itv]]][i]*/
-             }else{
-               cov[ioffset+ipos]*=covar[Tvar[TvarFind[itv]]][i];
-             }
            }
            iposold=ipos;
-           /* For products */
+           cov[ioffset+ipos]=cotvarv;
          }
-         /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */
-         /*   iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */
-         /*   cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */
-         /*   k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */
-         /*   cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */
-         /*   printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */
-         /* } */
-         /* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates *\/ */
-         /*   iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */
-         /*   /\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */
-         /*   cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; */
-         /* } */
          /* for products of time varying to be done */
          for (ii=1;ii<=nlstate+ndeath;ii++)
            for (j=1;j<=nlstate+ndeath;j++){
- Line 3984  double func( double *x)
+ Line 4089  double func( double *x)
            cov[2]=agexact;
            if(nagesqr==1)
              cov[3]= agexact*agexact;  /* Should be changed here */
-           for (kk=1; kk<=cptcovage;kk++) {
+           /* for (kk=1; kk<=cptcovage;kk++) { */
-             if(!FixedV[Tvar[Tage[kk]]])
+           /*   if(!FixedV[Tvar[Tage[kk]]]) */
-               cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
+           /*     cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /\* Tage[kk] gives the data-covariate associated with age *\/ */
-             else
+           /*   else */
-               cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact;
+           /*     cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\/  */
+           /* } */
+           for(ncovva=1, iposold=0; ncovva <= ncovta ; ncovva++){ /* Time varying  covariates with age including individual from products, product is computed dynamically */
+             itv=TvarAVVA[ncovva]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate, exploding product Vn*Vm into Vn and then Vm  */
+             ipos=TvarAVVAind[ncovva]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
+             if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */
+               cotvarv=cotvar[mw[mi][i]][TvarAVVA[ncovva]][i];  /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */
+             }else{ /* fixed covariate */
+               cotvarv=covar[itv][i];  /* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model */
+             }
+             if(ipos!=iposold){ /* Not a product or first of a product */
+               cotvarvold=cotvarv;
+             }else{ /* A second product */
+               cotvarv=cotvarv*cotvarvold;
+             }
+             iposold=ipos;
+             cov[ioffset+ipos]=cotvarv*agexact;
+             /* For products */
            }
            out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
 ,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
            savm=oldm;
- Line 4084  double func( double *x)
+ Line 4207  double func( double *x)
          }
          /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
          /*if(lli ==000.0)*/
-         /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
+         /* printf("num[i], i=%d, bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
          ipmx +=1;
          sw += weight[i];
          ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
- Line 4101  double func( double *x)
+ Line 4224  double func( double *x)
          cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];
        for(mi=1; mi<= wav[i]-1; mi++){
          for(k=1; k <= ncovv ; k++){
-           cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
+           cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
          }
          for (ii=1;ii<=nlstate+ndeath;ii++)
            for (j=1;j<=nlstate+ndeath;j++){
- Line 4148  double func( double *x)
+ Line 4271  double func( double *x)
            if(nagesqr==1)
              cov[3]= agexact*agexact;
            for (kk=1; kk<=cptcovage;kk++) {
-             cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+             if(!FixedV[Tvar[Tage[kk]]])
+               cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
+             else
+               cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
            }
            out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
 ,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- Line 4204  double func( double *x)
+ Line 4330  double func( double *x)
          ipmx +=1;
          sw += weight[i];
          ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
- /*      printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
+         /* printf("num[i]=%09ld, i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
        } /* end of wave */
      } /* end of individual */
    }else{  /* ml=5 no inter-extrapolation no jackson =0.8a */
- Line 4223  double func( double *x)
+ Line 4349  double func( double *x)
            if(nagesqr==1)
              cov[3]= agexact*agexact;
            for (kk=1; kk<=cptcovage;kk++) {
-             cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+             if(!FixedV[Tvar[Tage[kk]]])
+               cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
+             else
+               cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
            }
            out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- Line 4252  double func( double *x)
+ Line 4381  double func( double *x)
  double funcone( double *x)
  {
    /* Same as func but slower because of a lot of printf and if */
-   int i, ii, j, k, mi, d, kk, kf=0;
+   int i, ii, j, k, mi, d, kk, kv=0, kf=0;
    int ioffset=0;
    int ipos=0,iposold=0,ncovv=0;
+   double cotvarv, cotvarvold;
    double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
    double **out;
    double lli; /* Individual log likelihood */
- Line 4287  double funcone( double *x)
+ Line 4417  double funcone( double *x)
      /* Fixed */
      /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */
      /* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */
-     for (kf=1; kf<=ncovf;kf++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */
+     for (kf=1; kf<=ncovf;kf++){ /*  V2  +  V3  +  V4  Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */
        /* printf("Debug3 TvarFind[%d]=%d",kf, TvarFind[kf]); */
        /* printf(" Tvar[TvarFind[kf]]=%d", Tvar[TvarFind[kf]]); */
        /* printf(" i=%d covar[Tvar[TvarFind[kf]]][i]=%f\n",i,covar[Tvar[TvarFind[kf]]][i]); */
- Line 4310  double funcone( double *x)
+ Line 4440  double funcone( double *x)
           mw[mi][i] is real wave of the mi th effectve wave */
        /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i];
           s2=s[mw[mi+1][i]][i];
-          And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i]
+          And the iv th varying covariate in the DATA is the cotvar[mw[mi+1][i]][ncovcol+nqv+iv][i]
-          But if the variable is not in the model TTvar[iv] is the real variable effective in the model:
-          meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]
        */
      /* This part may be useless now because everythin should be in covar */
      /* for (k=1; k<=nqfveff;k++){ /\* Simple and product fixed Quantitative covariates without age* products *\/ */
- Line 4337  double funcone( double *x)
+ Line 4465  double funcone( double *x)
        /* We need the position of the time varying or product in the model */
        /* TvarVVind={2,5,5}, for V3 at position 2 and then the product V1*V3 is decomposed into V1 and V3 but at same position 5 */
        /* TvarVV gives the variable name */
-       for(ncovv=1, ipos=0; ncovv <= ncovvt ; ncovv++){ /* Varying  covariates (single and product but no age )*/
+       /* Other example V1 + V3 + V5 + age*V1  + age*V3 + age*V5 + V1*V3  + V3*V5  + V1*V5
-         itv=TvarVV[ncovv]; /*  TvarVV={3, 1, 3}  */
+       *      k=         1   2     3     4         5        6        7       8        9
-         ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] */
+       *  varying            1     2                                 3       4        5
+       *  ncovv              1     2                                3 4     5 6      7 8
+       * TvarVV[ncovv]      V3     5                                1 3     3 5      1 5
+       * TvarVVind           2     3                                7 7     8 8      9 9
+       * TvarFind[k]     1   0     0     0         0        0        0       0        0
+       */
+       /* Other model ncovcol=5 nqv=0 ntv=3 nqtv=0 nlstate=3
+        * V2 V3 V4 are fixed V6 V7 are timevarying so V8 and V5 are not in the model and product column will start at 9 Tvar[(v6*V2)6]=9
+         * FixedV[ncovcol+qv+ntv+nqtv]       V5
+         * 3           V1  V2     V3    V4   V5 V6     V7  V8 V3*V2 V7*V2  V6*V3 V7*V3 V6*V4 V7*V4
+         *             0   0      0      0    0  1      1   1  0, 0, 1,1,   1, 0, 1, 0, 1, 0, 1, 0}
+         * 3           0   0      0      0    0  1      1   1  0,     1      1    1      1    1}
+         * model=          V2  +  V3  +  V4  +  V6  +  V7  +  V6*V2  +  V7*V2  +  V6*V3  +  V7*V3  +  V6*V4  +  V7*V4
+         *                +age*V2 +age*V3 +age*V4 +age*V6 + age*V7
+         *                +age*V6*V2 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4
+         * model2=          V2  +  V3  +  V4  +  V6  +  V7  +  V3*V2  +  V7*V2  +  V6*V3  +  V7*V3  +  V6*V4  +  V7*V4
+         *                +age*V2 +age*V3 +age*V4 +age*V6 + age*V7
+         *                +age*V3*V2 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4
+         * model3=          V2  +  V3  +  V4  +  V6  +  V7  + age*V3*V2  +  V7*V2  +  V6*V3  +  V7*V3  +  V6*V4  +  V7*V4
+         *                +age*V2 +age*V3 +age*V4 +age*V6 + age*V7
+         *                +V3*V2 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4
+         * kmodel           1     2      3      4      5        6         7         8         9        10        11
+         *                  12       13      14      15       16
+         *                    17        18         19        20         21
+         * Tvar[kmodel]     2     3      4      6      7        9        10        11        12        13        14
+         *                   2       3        4       6        7
+         *                     9         11          12        13         14
+         * cptcovage=5+5 total of covariates with age
+         * Tage[cptcovage] age*V2=12      13      14      15       16
+         *1                   17            18         19        20         21 gives the position in model of covariates associated with age
+         *3 Tage[cptcovage] age*V3*V2=6
+         *3                age*V2=12         13      14      15       16
+         *3                age*V6*V3=18      19    20   21
+         * Tvar[Tage[cptcovage]]    Tvar[12]=2      3      4       6         Tvar[16]=7(age*V7)
+         *     Tvar[17]age*V6*V2=9      Tvar[18]age*V6*V3=11  age*V7*V3=12         age*V6*V4=13        Tvar[21]age*V7*V4=14
+         * 2   Tvar[17]age*V3*V2=9      Tvar[18]age*V6*V3=11  age*V7*V3=12         age*V6*V4=13        Tvar[21]age*V7*V4=14
+         * 3 Tvar[Tage[cptcovage]]    Tvar[6]=9      Tvar[12]=2      3     4       6         Tvar[16]=7(age*V7)
+         * 3     Tvar[18]age*V6*V3=11  age*V7*V3=12         age*V6*V4=13        Tvar[21]age*V7*V4=14
+         * 3 Tage[cptcovage] age*V3*V2=6   age*V2=12 age*V3 13    14      15       16
+         *                    age*V6*V3=18         19        20         21 gives the position in model of covariates associated with age
+         * 3   Tvar[17]age*V3*V2=9      Tvar[18]age*V6*V3=11  age*V7*V3=12         age*V6*V4=13        Tvar[21]age*V7*V4=14
+         * Tvar=                {2, 3, 4, 6, 7,
+         *                       9, 10, 11, 12, 13, 14,
+         *              Tvar[12]=2, 3, 4, 6, 7,
+         *              Tvar[17]=9, 11, 12, 13, 14}
+         * Typevar[1]@21 = {0, 0, 0, 0, 0,
+         *                  2, 2, 2, 2, 2, 2,
+         * 3                3, 2, 2, 2, 2, 2,
+         *                  1, 1, 1, 1, 1,
+         *                  3, 3, 3, 3, 3}
+         * 3                 2, 3, 3, 3, 3}
+         * p Tposprod[1]@21 {0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 1, 3, 4, 5, 6} Id of the prod at position k in the model
+         * p Tprod[1]@21 {6, 7, 8, 9, 10, 11, 0 <repeats 15 times>}
+         * 3 Tposprod[1]@21 {0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 1, 3, 4, 5, 6}
+         * 3 Tprod[1]@21 {17, 7, 8, 9, 10, 11, 0 <repeats 15 times>}
+         * cptcovprod=11 (6+5)
+         * FixedV[Tvar[Tage[cptcovage]]]]  FixedV[2]=0      FixedV[3]=0      0      1          (age*V7)Tvar[16]=1 FixedV[absolute] not [kmodel]
+         *   FixedV[Tvar[17]=FixedV[age*V6*V2]=FixedV[9]=1        1         1          1         1
+         * 3 FixedV[Tvar[17]=FixedV[age*V3*V2]=FixedV[9]=0        [11]=1         1          1         1
+         * FixedV[]          V1=0     V2=0   V3=0  v4=0    V5=0  V6=1    V7=1 v8=1  OK then model dependent
+         *                   9=1  [V7*V2]=[10]=1 11=1  12=1  13=1  14=1
+         * 3                 9=0  [V7*V2]=[10]=1 11=1  12=1  13=1  14=1
+         * cptcovdageprod=5  for gnuplot printing
+         * cptcovprodvage=6
+         * ncova=15           1        2       3       4       5
+         *                      6 7        8 9      10 11        12 13     14 15
+         * TvarA              2        3       4       6       7
+         *                      6 2        6 7       7 3          6 4       7 4
+         * TvaAind             12 12      13 13     14 14      15 15       16 16
+         * ncovf            1     2      3
+         *                                    V6       V7      V6*V2     V7*V2     V6*V3     V7*V3     V6*V4     V7*V4
+         * ncovvt=14                            1      2        3 4       5 6       7 8       9 10     11 12     13 14
+         * TvarVV[1]@14 = itv               {V6=6,     7, V6*V2=6, 2,     7, 2,     6, 3,     7, 3,     6, 4,     7, 4}
+         * TvarVVind[1]@14=                    {4,     5,       6, 6,     7, 7,     8, 8,      9, 9,   10, 10,   11, 11}
+         * 3 ncovvt=12                        V6       V7      V7*V2     V6*V3     V7*V3     V6*V4     V7*V4
+         * 3 TvarVV[1]@12 = itv                {6,     7, V7*V2=7, 2,     6, 3,     7, 3,     6, 4,     7, 4}
+         * 3                                    1      2        3  4      5  6      7  8      9 10     11 12
+         * TvarVVind[1]@12=         {V6 is in k=4,     5,  7,(4isV2)=7,   8, 8,      9, 9,   10,10,    11,11}TvarVVind[12]=k=11
+         * TvarV              6, 7, 9, 10, 11, 12, 13, 14
+         * 3 cptcovprodvage=6
+         * 3 ncovta=15    +age*V3*V2+age*V2+agev3+ageV4 +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4
+         * 3 TvarAVVA[1]@15= itva 3 2    2      3    4        6       7        6 3         7 3         6 4         7 4
+         * 3 ncovta             1 2      3      4    5        6       7        8 9       10 11       12 13        14 15
+         * TvarAVVAind[1]@15= V3 is in k=2 1 1  2    3        4       5        4,2         5,2,      4,3           5 3}TvarVVAind[]
+         * TvarAVVAind[1]@15= V3 is in k=6 6 12  13   14      15      16       18 18       19,19,     20,20        21,21}TvarVVAind[]
+         * 3 ncovvta=10     +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4
+         * 3 we want to compute =cotvar[mw[mi][i]][TvarVVA[ncovva]][i] at position TvarVVAind[ncovva]
+         * 3 TvarVVA[1]@10= itva   6       7        6 3         7 3         6 4         7 4
+         * 3 ncovva                1       2        3 4         5 6         7 8         9 10
+         * TvarVVAind[1]@10= V6 is in k=4  5        8,8         9, 9,      10,10        11 11}TvarVVAind[]
+         * TvarVVAind[1]@10=       15       16     18,18        19,19,      20,20        21 21}TvarVVAind[]
+         * TvarVA              V3*V2=6 6 , 1, 2, 11, 12, 13, 14
+         * TvarFind[1]@14= {1,    2,     3,     0 <repeats 12 times>}
+         * Tvar[1]@21=     {2,    3,     4,    6,      7,      9,      10,        11,       12,      13,       14,
+         *                   2, 3, 4, 6, 7,
+         *                     6, 8, 9, 10, 11}
+         * TvarFind[itv]                        0      0       0
+         * FixedV[itv]                          1      1       1  0      1 0       1 0       1 0       0
+         * Tvar[TvarFind[ncovf]]=[1]=2 [2]=3 [4]=4
+         * Tvar[TvarFind[itv]]                [0]=?      ?ncovv 1 à ncovvt]
+         *   Not a fixed cotvar[mw][itv][i]     6       7      6  2      7, 2,     6, 3,     7, 3,     6, 4,     7, 4}
+         *   fixed covar[itv]                  [6]     [7]    [6][2]
+         */
+       for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /*  V6       V7      V7*V2     V6*V3     V7*V3     V6*V4     V7*V4 Time varying  covariates (single and extended product but no age) including individual from products, product is computed dynamically */
+         itv=TvarVV[ncovv]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate, or fixed covariate of a varying product after exploding product Vn*Vm into Vn and then Vm  */
+         ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
+         /* if(TvarFind[itv]==0){ /\* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv *\/ */
+         if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */
+           cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i];  /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */
+         }else{ /* fixed covariate */
+           /* cotvarv=covar[Tvar[TvarFind[itv]]][i];  /\* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model *\/ */
+           cotvarv=covar[itv][i];  /* Good: In V6*V3, 3 is fixed at position of the data */
+         }
          if(ipos!=iposold){ /* Not a product or first of a product */
-           /* TvarFind={1,0,0,0}  */
+           cotvarvold=cotvarv;
-           if(TvarFind[itv]==0){
+         }else{ /* A second product */
-             cov[ioffset+ipos]= cotvar[mw[mi][i]][ncovv][i];  /* Should be covar if fixed covar[Tvar[TvarFind[itv]]][i]*/
+           cotvarv=cotvarv*cotvarvold;
-           }else{
-             cov[ioffset+ipos]=covar[Tvar[TvarFind[itv]]][i];
-           }
-         }else{
-           if(TvarFind[itv]==0){
-             cov[ioffset+ipos]*= cotvar[mw[mi][i]][ncovv][i];  /* Should be covar if fixed covar[Tvar[TvarFind[itv]]][i]*/
-           }else{
-             cov[ioffset+ipos]*=covar[Tvar[TvarFind[itv]]][i];
-           }
          }
          iposold=ipos;
+         cov[ioffset+ipos]=cotvarv;
          /* For products */
        }
        /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates single *\/ */
- Line 4393  double funcone( double *x)
+ Line 4626  double funcone( double *x)
          cov[2]=agexact;
          if(nagesqr==1)
            cov[3]= agexact*agexact;
-         for (kk=1; kk<=cptcovage;kk++) {
+         for(ncovva=1, iposold=0; ncovva <= ncovta ; ncovva++){ /* Time varying  covariates with age including individual from products, product is computed dynamically */
-           if(!FixedV[Tvar[Tage[kk]]])
+           itv=TvarAVVA[ncovva]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate, exploding product Vn*Vm into Vn and then Vm  */
-             cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
+           ipos=TvarAVVAind[ncovva]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
-           else
+           /* if(TvarFind[itv]==0){ /\* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv *\/ */
-             cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact;
+           if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */
+             /* printf("DEBUG  ncovva=%d, Varying TvarAVVA[ncovva]=%d\n", ncovva, TvarAVVA[ncovva]); */
+             cotvarv=cotvar[mw[mi][i]][TvarAVVA[ncovva]][i];  /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */
+           }else{ /* fixed covariate */
+             /* cotvarv=covar[Tvar[TvarFind[itv]]][i];  /\* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model *\/ */
+             /* printf("DEBUG ncovva=%d, Fixed TvarAVVA[ncovva]=%d\n", ncovva, TvarAVVA[ncovva]); */
+             cotvarv=covar[itv][i];  /* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model */
+           }
+           if(ipos!=iposold){ /* Not a product or first of a product */
+             cotvarvold=cotvarv;
+           }else{ /* A second product */
+             /* printf("DEBUG * \n"); */
+             cotvarv=cotvarv*cotvarvold;
+           }
+           iposold=ipos;
+           /* printf("DEBUG Product cov[ioffset+ipos=%d] \n",ioffset+ipos); */
+           cov[ioffset+ipos]=cotvarv*agexact;
+           /* For products */
          }
          /* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */
          /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
          out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
- Line 4452  double funcone( double *x)
+ Line 4703  double funcone( double *x)
        ipmx +=1;
        sw += weight[i];
        ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
-       /* printf("Funcone i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
+       /* Printing covariates values for each contribution for checking */
+       /* printf("num[i]=%09ld, i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
        if(globpr){
          fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
   %11.6f %11.6f %11.6f ", \
                  num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
 *weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2]));
-  /*     printf("%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ */
+         /*      printf("%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ */
-  /* %11.6f %11.6f %11.6f ", \ */
+         /* %11.6f %11.6f %11.6f ", \ */
-  /*             num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, */
+         /*              num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, */
-  /*             2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
+         /*              2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
          for(k=1,llt=0.,l=0.; k<=nlstate; k++){
            llt +=ll[k]*gipmx/gsw;
            fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
            /* printf(" %10.6f",-ll[k]*gipmx/gsw); */
          }
-         fprintf(ficresilk," %10.6f\n", -llt);
+         fprintf(ficresilk," %10.6f ", -llt);
          /* printf(" %10.6f\n", -llt); */
-       }
+         /* if(debugILK){ /\* debugILK is set by a #d in a comment line *\/ */
+         /* fprintf(ficresilk,"%09ld ", num[i]); */ /* not necessary */
+         for (kf=1; kf<=ncovf;kf++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */
+           fprintf(ficresilk," %g",covar[Tvar[TvarFind[kf]]][i]);
+         }
+         for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Varying  covariates (single and product but no age) including individual from products */
+           ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
+           if(ipos!=iposold){ /* Not a product or first of a product */
+             fprintf(ficresilk," %g",cov[ioffset+ipos]);
+             /* printf(" %g",cov[ioffset+ipos]); */
+           }else{
+             fprintf(ficresilk,"*");
+             /* printf("*"); */
+           }
+           iposold=ipos;
+         }
+         /* for (kk=1; kk<=cptcovage;kk++) { */
+         /*   if(!FixedV[Tvar[Tage[kk]]]){ */
+         /*     fprintf(ficresilk," %g*age",covar[Tvar[Tage[kk]]][i]); */
+         /*     /\* printf(" %g*age",covar[Tvar[Tage[kk]]][i]); *\/ */
+         /*   }else{ */
+         /*     fprintf(ficresilk," %g*age",cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]);/\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\/  */
+         /*     /\* printf(" %g*age",cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]);/\\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\\/  *\/ */
+         /*   } */
+         /* } */
+         for(ncovva=1, iposold=0; ncovva <= ncovta ; ncovva++){ /* Time varying  covariates with age including individual from products, product is computed dynamically */
+           itv=TvarAVVA[ncovva]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate, exploding product Vn*Vm into Vn and then Vm  */
+           ipos=TvarAVVAind[ncovva]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
+           /* if(TvarFind[itv]==0){ /\* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv *\/ */
+           if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */
+             /* printf("DEBUG  ncovva=%d, Varying TvarAVVA[ncovva]=%d\n", ncovva, TvarAVVA[ncovva]); */
+             cotvarv=cotvar[mw[mi][i]][TvarAVVA[ncovva]][i];  /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */
+           }else{ /* fixed covariate */
+             /* cotvarv=covar[Tvar[TvarFind[itv]]][i];  /\* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model *\/ */
+             /* printf("DEBUG ncovva=%d, Fixed TvarAVVA[ncovva]=%d\n", ncovva, TvarAVVA[ncovva]); */
+             cotvarv=covar[itv][i];  /* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model */
+           }
+           if(ipos!=iposold){ /* Not a product or first of a product */
+             cotvarvold=cotvarv;
+           }else{ /* A second product */
+             /* printf("DEBUG * \n"); */
+             cotvarv=cotvarv*cotvarvold;
+           }
+           cotvarv=cotvarv*agexact;
+           fprintf(ficresilk," %g*age",cotvarv);
+           iposold=ipos;
+           /* printf("DEBUG Product cov[ioffset+ipos=%d] \n",ioffset+ipos); */
+           cov[ioffset+ipos]=cotvarv;
+           /* For products */
+         }
+         /* printf("\n"); */
+         /* } /\*  End debugILK *\/ */
+         fprintf(ficresilk,"\n");
+       } /* End if globpr */
      } /* end of wave */
    } /* end of individual */
    for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
- Line 4479  double funcone( double *x)
+ Line 4784  double funcone( double *x)
      gipmx=ipmx;
      gsw=sw;
    }
- return -l;
+   return -l;
  }
- Line 4490  void likelione(FILE *ficres,double p[],
+ Line 4795  void likelione(FILE *ficres,double p[],
       the selection of individuals/waves and
       to check the exact contribution to the likelihood.
       Plotting could be done.
-    */
+   */
-   int k;
+   void pstamp(FILE *ficres);
+   int k, kf, kk, kvar, ncovv, iposold, ipos;
    if(*globpri !=0){ /* Just counts and sums, no printings */
      strcpy(fileresilk,"ILK_");
- Line 4500  void likelione(FILE *ficres,double p[],
+ Line 4806  void likelione(FILE *ficres,double p[],
        printf("Problem with resultfile: %s\n", fileresilk);
        fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);
      }
+     pstamp(ficresilk);fprintf(ficresilk,"# model=1+age+%s\n",model);
      fprintf(ficresilk, "#individual(line's_record) count ageb ageend s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
      fprintf(ficresilk, "#num_i ageb agend i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav ");
      /*  i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */
      for(k=1; k<=nlstate; k++)
        fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);
-     fprintf(ficresilk," -2*gipw/gsw*weight*ll(total)\n");
+     fprintf(ficresilk," -2*gipw/gsw*weight*ll(total) ");
-   }
+     /* if(debugILK){ /\* debugILK is set by a #d in a comment line *\/ */
+       for(kf=1;kf <= ncovf; kf++){
+         fprintf(ficresilk,"V%d",Tvar[TvarFind[kf]]);
+         /* printf("V%d",Tvar[TvarFind[kf]]); */
+       }
+       for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){
+         ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate */
+         if(ipos!=iposold){ /* Not a product or first of a product */
+           /* printf(" %d",ipos); */
+           fprintf(ficresilk," V%d",TvarVV[ncovv]);
+         }else{
+           /* printf("*"); */
+           fprintf(ficresilk,"*");
+         }
+         iposold=ipos;
+       }
+       for (kk=1; kk<=cptcovage;kk++) {
+         if(!FixedV[Tvar[Tage[kk]]]){
+           /* printf(" %d*age(Fixed)",Tvar[Tage[kk]]); */
+           fprintf(ficresilk," %d*age(Fixed)",Tvar[Tage[kk]]);
+         }else{
+           fprintf(ficresilk," %d*age(Varying)",Tvar[Tage[kk]]);/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
+           /* printf(" %d*age(Varying)",Tvar[Tage[kk]]);/\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\/  */
+         }
+       }
+     /* } /\* End if debugILK *\/ */
+     /* printf("\n"); */
+     fprintf(ficresilk,"\n");
+   } /* End glogpri */
    *fretone=(*func)(p);
    if(*globpri !=0){
- Line 4518  void likelione(FILE *ficres,double p[],
+ Line 4854  void likelione(FILE *ficres,double p[],
      fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
      fprintf(fichtm,"\n<br>Equation of the model: <b>model=1+age+%s</b><br>\n",model);
-     for (k=1; k<= nlstate ; k++) {
-       fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
- <img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k);
-     }
      fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \
- <img src=\"%s-ori.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
+ <img src=\"%s-ori.png\">\n",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
      fprintf(fichtm,"<br>- and by state of destination <a href=\"%s-dest.png\">%s-dest.png</a><br> \
- <img src=\"%s-dest.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
+ <img src=\"%s-dest.png\">\n",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
+     for (k=1; k<= nlstate ; k++) {
+       fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br>\n \
+ <img src=\"%s-p%dj.png\">\n",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k);
+       for(kf=1; kf <= ncovf; kf++){ /* For each simple dummy covariate of the model */
+          kvar=Tvar[TvarFind[kf]];  /* variable */
+          fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j with colored covariate V%d. Same dot size of all points but with a different color for transitions with dummy variable V%d=1 at beginning of transition (keeping former color for V%d=0): ",k,k,Tvar[TvarFind[kf]],Tvar[TvarFind[kf]],Tvar[TvarFind[kf]]);
+          fprintf(fichtm,"<a href=\"%s-p%dj-%d.png\">%s-p%dj-%d.png</a><br>",subdirf2(optionfilefiname,"ILK_"),k,kvar,subdirf2(optionfilefiname,"ILK_"),k,kvar);
+          fprintf(fichtm,"<img src=\"%s-p%dj-%d.png\">",subdirf2(optionfilefiname,"ILK_"),k,Tvar[TvarFind[kf]]);
+       }
+       for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Loop on the time varying extended covariates (with extension of Vn*Vm */
+         ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate */
+         kvar=TvarVV[ncovv]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate */
+         /* printf("DebugILK fichtm ncovv=%d, kvar=TvarVV[ncovv]=V%d, ipos=TvarVVind[ncovv]=%d, Dummy[ipos]=%d, Typevar[ipos]=%d\n", ncovv,kvar,ipos,Dummy[ipos],Typevar[ipos]); */
+         if(ipos!=iposold){ /* Not a product or first of a product */
+           /* fprintf(ficresilk," V%d",TvarVV[ncovv]); */
+           /* printf(" DebugILK fichtm ipos=%d != iposold=%d\n", ipos, iposold); */
+           if(Dummy[ipos]==0 && Typevar[ipos]==0){ /* Only if dummy time varying: Dummy(0, 1=quant singor prod without age,2 dummy*age, 3quant*age) Typevar (0 single, 1=*age,2=Vn*vm)  */
+             fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j with colored time varying dummy covariate V%d. Same dot size of all points but with a different color for transitions with dummy variable V%d=1 at beginning of transition (keeping former color for V%d=0): <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
+ <img src=\"%s-p%dj-%d.png\">",k,k,kvar,kvar,kvar,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,kvar);
+           } /* End only for dummies time varying (single?) */
+         }else{ /* Useless product */
+           /* printf("*"); */
+           /* fprintf(ficresilk,"*"); */
+         }
+         iposold=ipos;
+       } /* For each time varying covariate */
+     } /* End loop on states */
+ /*     if(debugILK){ */
+ /*       for(kf=1; kf <= ncovf; kf++){ /\* For each simple dummy covariate of the model *\/ */
+ /*      /\* kvar=Tvar[TvarFind[kf]]; *\/ /\* variable *\/ */
+ /*      for (k=1; k<= nlstate ; k++) { */
+ /*        fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j with colored covariate V%. Same dot size of all points but with a different color for transitions with dummy variable V%d=1 at beginning of transition (keeping former color for V%d=0): <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \ */
+ /* <img src=\"%s-p%dj-%d.png\">",k,k,Tvar[TvarFind[kf]],Tvar[TvarFind[kf]],Tvar[TvarFind[kf]],subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,Tvar[TvarFind[kf]]); */
+ /*      } */
+ /*       } */
+ /*       for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /\* Loop on the time varying extended covariates (with extension of Vn*Vm *\/ */
+ /*      ipos=TvarVVind[ncovv]; /\* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate *\/ */
+ /*      kvar=TvarVV[ncovv]; /\*  TvarVV={3, 1, 3} gives the name of each varying covariate *\/ */
+ /*      /\* printf("DebugILK fichtm ncovv=%d, kvar=TvarVV[ncovv]=V%d, ipos=TvarVVind[ncovv]=%d, Dummy[ipos]=%d, Typevar[ipos]=%d\n", ncovv,kvar,ipos,Dummy[ipos],Typevar[ipos]); *\/ */
+ /*      if(ipos!=iposold){ /\* Not a product or first of a product *\/ */
+ /*        /\* fprintf(ficresilk," V%d",TvarVV[ncovv]); *\/ */
+ /*        /\* printf(" DebugILK fichtm ipos=%d != iposold=%d\n", ipos, iposold); *\/ */
+ /*        if(Dummy[ipos]==0 && Typevar[ipos]==0){ /\* Only if dummy time varying: Dummy(0, 1=quant singor prod without age,2 dummy*age, 3quant*age) Typevar (0 single, 1=*age,2=Vn*vm)  *\/ */
+ /*          for (k=1; k<= nlstate ; k++) { */
+ /*            fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Same dot size of all points but with a different color for transitions with dummy variable V%d=1 at beginning of transition (keeping former color for V%d=0): <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \ */
+ /* <img src=\"%s-p%dj-%d.png\">",k,k,kvar,kvar,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,kvar); */
+ /*          } /\* End state *\/ */
+ /*        } /\* End only for dummies time varying (single?) *\/ */
+ /*      }else{ /\* Useless product *\/ */
+ /*        /\* printf("*"); *\/ */
+ /*        /\* fprintf(ficresilk,"*"); *\/  */
+ /*      } */
+ /*      iposold=ipos; */
+ /*       } /\* For each time varying covariate *\/ */
+ /*     }/\* End debugILK *\/ */
      fflush(fichtm);
-   }
+   }/* End globpri */
    return;
  }
- Line 5249  Title=%s <br>Datafile=%s Firstpass=%d La
+ Line 5638  Title=%s <br>Datafile=%s Firstpass=%d La
                if(anyvaryingduminmodel==1){ /* Some are varying covariates */
                  for (z1=1; z1<=cptcoveff; z1++) {
                    if( Fixed[Tmodelind[z1]]==1){
-                     iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
+                     /* iv= Tvar[Tmodelind[z1]]-ncovcol-nqv; /\* Good *\/ */
+                     iv= Tvar[Tmodelind[z1]]; /* Good *//* because cotvar starts now at first at ncovcol+nqv+ntv */
                      if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[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. */
- Line 5330  Title=%s <br>Datafile=%s Firstpass=%d La
+ Line 5720  Title=%s <br>Datafile=%s Firstpass=%d La
          fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");
          fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");
          fprintf(ficlog, "\n#********** Variable ");
-         for (z1=1; z1<=cptcovs; z1++){
+         for (z1=1; z1<=cptcoveff; z1++){
            if(!FixedV[Tvaraff[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]])]);
- Line 5786  void prevalence(double ***probs, double
+ Line 6176  void prevalence(double ***probs, double
          /* Tvar[Tmodelind[z1]] is the n of Vn; n-ncovcol-nqv is the first time varying covariate or iv */
          for (z1=1; z1<=cptcoveff; z1++){
            if( Fixed[Tmodelind[z1]]==1){
-             iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
+             iv= Tvar[Tmodelind[z1]];/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
              if (cotvar[m][iv][i]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality */
                bool=0;
            }else if( Fixed[Tmodelind[z1]]== 0)  /* fixed */
- Line 6095  void  concatwav(int wav[], int **dh, int
+ Line 6485  void  concatwav(int wav[], int **dh, int
     /* Loop on covariates without age and products and no quantitative variable */
     for (k=1; k<=cptcovt; k++) { /* cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */
       for (j=-1; (j < maxncov); j++) Ndum[j]=0;
-      printf("Testing k=%d, cptcovt=%d\n",k, cptcovt);
+      /* printf("Testing k=%d, cptcovt=%d\n",k, cptcovt); */
-      if(Dummy[k]==0 && Typevar[k] !=1 && Typevar[k] != 2){ /* Dummy covariate and not age product nor fixed product */
+      if(Dummy[k]==0 && Typevar[k] !=1 && Typevar[k] != 3  && Typevar[k] != 2){ /* Dummy covariate and not age product nor fixed product */
         switch(Fixed[k]) {
         case 0: /* Testing on fixed dummy covariate, simple or product of fixed */
           modmaxcovj=0;
- Line 6193  void  concatwav(int wav[], int **dh, int
+ Line 6583  void  concatwav(int wav[], int **dh, int
           break;
         } /* end switch */
       } /* end dummy test */
-      if(Dummy[k]==1 && Typevar[k] !=1){ /* Quantitative covariate and not age product */
+      if(Dummy[k]==1 && Typevar[k] !=1 && Typevar[k] !=3 && Fixed ==0){ /* Fixed Quantitative covariate and not age product */
         for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the  modality of this covariate Vj*/
           if(Tvar[k]<=0 || Tvar[k]>=NCOVMAX){
             printf("Error k=%d \n",k);
- Line 6659  void  concatwav(int wav[], int **dh, int
+ Line 7049  void  concatwav(int wav[], int **dh, int
     /*   fprintf(ficresprobmorprev," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */
     /* } */
     for (j=1; j<= cptcovs; j++){ /* For each selected (single) quantitative value */ /* To be done*/
+      /* fprintf(ficresprobmorprev," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); */
       fprintf(ficresprobmorprev," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
     }
     /* for(j=1;j<=cptcoveff;j++)  */
- Line 7293  To be simple, these graphs help to under
+ Line 7684  To be simple, these graphs help to under
     for(nres=1;nres <=nresult; nres++){ /* For each resultline */
     for(j1=1; j1<=tj;j1++){ /* For any combination of dummy covariates, fixed and varying */
-      printf("Varprob  TKresult[nres]=%d j1=%d, nres=%d, cptcovn=%d, cptcoveff=%d tj=%d cptcovs=%d\n",  TKresult[nres], j1, nres, cptcovn, cptcoveff, tj, cptcovs);
+      /* printf("Varprob  TKresult[nres]=%d j1=%d, nres=%d, cptcovn=%d, cptcoveff=%d tj=%d cptcovs=%d\n",  TKresult[nres], j1, nres, cptcovn, cptcoveff, tj, cptcovs); */
       if(tj != 1 && TKresult[nres]!= j1)
         continue;
- Line 7309  To be simple, these graphs help to under
+ Line 7700  To be simple, these graphs help to under
         /* Including quantitative variables of the resultline to be done */
         for (z1=1; z1<=cptcovs; z1++){ /* Loop on each variable of this resultline  */
-          printf("Varprob modelresult[%d][%d]=%d model=1+age+%s \n",nres, z1, modelresult[nres][z1], model);
+          /* printf("Varprob modelresult[%d][%d]=%d model=1+age+%s \n",nres, z1, modelresult[nres][z1], model); */
           fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=1+age+%s \n",nres, z1, modelresult[nres][z1], model);
           /* fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=1+age+%s resultline[%d]=%s \n",nres, z1, modelresult[nres][z1], model, nres, resultline[nres]); */
           if(Dummy[modelresult[nres][z1]]==0){/* Dummy variable of the variable in position modelresult in the model corresponding to z1 in resultline  */
- Line 7375  To be simple, these graphs help to under
+ Line 7766  To be simple, these graphs help to under
           cov[3]= age*age;
         /* New code end of combination but for each resultline */
         for(k1=1;k1<=cptcovt;k1++){ /* loop on model equation (including products) */
-          if(Typevar[k1]==1){ /* A product with age */
+          if(Typevar[k1]==1 || Typevar[k1] ==3){ /* A product with age */
             cov[2+nagesqr+k1]=precov[nres][k1]*cov[2];
           }else{
             cov[2+nagesqr+k1]=precov[nres][k1];
- Line 7968  void printinggnuplot(char fileresu[], ch
+ Line 8359  void printinggnuplot(char fileresu[], ch
    char dirfileres[132],optfileres[132];
    char gplotcondition[132], gplotlabel[132];
-   int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,ij=0, ijp=0, l=0;
+   int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,kf=0,kvar=0,kk=0,ipos=0,iposold=0,ij=0, ijp=0, l=0;
    int lv=0, vlv=0, kl=0;
    int ng=0;
    int vpopbased;
- Line 7994  void printinggnuplot(char fileresu[], ch
+ Line 8385  void printinggnuplot(char fileresu[], ch
    fprintf(ficgp,"yoff=(%d > 2? 0:1);\n",nlstate);
    fprintf(ficgp,"\n#Peripheral arrows\nset for [i=1:%d] for [j=1:%d] arrow i*10+j from cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.95*(cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0) - cos(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta2:0)), -0.95*(sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) - sin(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d))+( i!=j?(i-j)/abs(i-j)*delta2:0)) ls (i < j? 1:2)\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
-   fprintf(ficgp,"\n#Centripete arrows (turning in other direction (1-i) instead of (i-1)) \nset for [i=1:%d] arrow (%d+1)*10+i from cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.80*(cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0)  ), -0.80*(sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) + yoff ) ls 4\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
+   fprintf(ficgp,"\n#Centripete arrows (turning in other direction (1-i) instead of (i-1)) \nset for [i=1:%d] for [j=1:%d] arrow (%d+1)*10+i from cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.80*(cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0)  ), -0.80*(sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) + yoff ) ls 4\n",nlstate, nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
    fprintf(ficgp,"\n#show arrow\nunset label\n");
    fprintf(ficgp,"\n#States labels, starting from 2 (2-i) instead of (1-i), was (i-1)\nset for [i=1:%d] label i sprintf(\"State %%d\",i) center at cos(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)), yoff+sin(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)) font \"helvetica, 16\" tc rgbcolor \"blue\"\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
    fprintf(ficgp,"\nset label %d+1 sprintf(\"State %%d\",%d+1) center at 0.,0.  font \"helvetica, 16\" tc rgbcolor \"red\"\n",nlstate,nlstate);
- Line 8031  void printinggnuplot(char fileresu[], ch
+ Line 8422  void printinggnuplot(char fileresu[], ch
    fprintf(ficgp,"\nset out;unset log\n");
    /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
+   /* Plot the probability implied in the likelihood by covariate value */
+   fprintf(ficgp,"\nset ter pngcairo size 640, 480");
+   /* if(debugILK==1){ */
+   for(kf=1; kf <= ncovf; kf++){ /* For each simple dummy covariate of the model */
+     kvar=Tvar[TvarFind[kf]]; /* variable name */
+     /* k=18+Tvar[TvarFind[kf]];/\*offset because there are 18 columns in the ILK_ file but could be placed else where *\/ */
+     /* k=18+kf;/\*offset because there are 18 columns in the ILK_ file *\/ */
+     k=19+kf;/*offset because there are 19 columns in the ILK_ file */
+     for (i=1; i<= nlstate ; i ++) {
+       fprintf(ficgp,"\nset out \"%s-p%dj-%d.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i,kvar);
+       fprintf(ficgp,"unset log;\n# For each simple dummy covariate of the model \n plot  \"%s\"",subdirf(fileresilk));
+       if(gnuplotversion >=5.2){ /* Former gnuplot versions do not have variable pointsize!! */
+         fprintf(ficgp,"  u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? 7 : 9):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt variable ps 0.4 lc variable \\\n",i,1,k,k,i,1,kvar);
+         for (j=2; j<= nlstate+ndeath ; j ++) {
+           fprintf(ficgp,",\\\n \"\" u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? 7 : 9):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt variable ps 0.4 lc variable ",i,j,k,k,i,j,kvar);
+         }
+       }else{
+         fprintf(ficgp,"  u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt 7 ps 0.4 lc variable \\\n",i,1,k,i,1,kvar);
+         for (j=2; j<= nlstate+ndeath ; j ++) {
+           fprintf(ficgp,",\\\n \"\" u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt 7 ps 0.4 lc variable ",i,j,k,i,j,kvar);
+         }
+       }
+       fprintf(ficgp,";\nset out; unset ylabel;\n");
+     }
+   } /* End of each covariate dummy */
+   for(ncovv=1, iposold=0, kk=0; ncovv <= ncovvt ; ncovv++){
+     /* Other example        V1 + V3 + V5 + age*V1  + age*V3 + age*V5 + V1*V3  + V3*V5  + V1*V5
+      *     kmodel       =     1   2     3     4         5        6        7       8        9
+      *  varying                   1     2                                 3       4        5
+      *  ncovv                     1     2                                3 4     5 6      7 8
+      * TvarVV[ncovv]             V3     5                                1 3     3 5      1 5
+      * TvarVVind[ncovv]=kmodel    2     3                                7 7     8 8      9 9
+      * TvarFind[kmodel]       1   0     0     0         0        0        0       0        0
+      * kdata     ncovcol=[V1 V2] nqv=0 ntv=[V3 V4] nqtv=V5
+      * Dummy[kmodel]          0   0     1     2         2        3        1       1        1
+      */
+     ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate */
+     kvar=TvarVV[ncovv]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate */
+     /* printf("DebugILK ficgp ncovv=%d, kvar=TvarVV[ncovv]=%d, ipos=TvarVVind[ncovv]=%d, Dummy[ipos]=%d, Typevar[ipos]=%d\n", ncovv,kvar,ipos,Dummy[ipos],Typevar[ipos]); */
+     if(ipos!=iposold){ /* Not a product or first of a product */
+       /* printf(" %d",ipos); */
+       /* fprintf(ficresilk," V%d",TvarVV[ncovv]); */
+       /* printf(" DebugILK ficgp suite ipos=%d != iposold=%d\n", ipos, iposold); */
+       kk++; /* Position of the ncovv column in ILK_ */
+       k=18+ncovf+kk; /*offset because there are 18 columns in the ILK_ file plus ncovf fixed covariate */
+       if(Dummy[ipos]==0 && Typevar[ipos]==0){ /* Only if dummy time varying: Dummy(0, 1=quant singor prod without age,2 dummy*age, 3quant*age) Typevar (0 single, 1=*age,2=Vn*vm)  */
+         for (i=1; i<= nlstate ; i ++) {
+           fprintf(ficgp,"\nset out \"%s-p%dj-%d.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i,kvar);
+           fprintf(ficgp,"unset log;\n# For each simple dummy covariate of the model \n plot  \"%s\"",subdirf(fileresilk));
+             /* printf("Before DebugILK gnuplotversion=%g >=5.2\n",gnuplotversion); */
+           if(gnuplotversion >=5.2){ /* Former gnuplot versions do not have variable pointsize!! */
+             /* printf("DebugILK gnuplotversion=%g >=5.2\n",gnuplotversion); */
+             fprintf(ficgp,"  u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? 7 : 9):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt variable ps 0.4 lc variable \\\n",i,1,k,k,i,1,kvar);
+             for (j=2; j<= nlstate+ndeath ; j ++) {
+               fprintf(ficgp,",\\\n \"\" u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? 7 : 9):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt variable ps 0.4 lc variable ",i,j,k,k,i,j,kvar);
+             }
+           }else{
+             /* printf("DebugILK gnuplotversion=%g <5.2\n",gnuplotversion); */
+             fprintf(ficgp,"  u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt 7 ps 0.4 lc variable \\\n",i,1,k,i,1,kvar);
+             for (j=2; j<= nlstate+ndeath ; j ++) {
+               fprintf(ficgp,",\\\n \"\" u  2:($5 == %d && $6==%d ? $10 : 1/0):($%d==0 ? $6 : $6+4) t \"p%d%d V%d\" with points pt 7 ps 0.4 lc variable ",i,j,k,i,j,kvar);
+             }
+           }
+           fprintf(ficgp,";\nset out; unset ylabel;\n");
+         }
+       }/* End if dummy varying */
+     }else{ /*Product */
+       /* printf("*"); */
+       /* fprintf(ficresilk,"*"); */
+     }
+     iposold=ipos;
+   } /* For each time varying covariate */
+   /* } /\* debugILK==1 *\/ */
+   /* unset log; plot  "rrtest1_sorted_4/ILK_rrtest1_sorted_4.txt" u  2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with points lc variable */
+   /* fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); */
+   /* fprintf(ficgp,"\nreplot  \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); */
+   fprintf(ficgp,"\nset out;unset log\n");
+   /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
    strcpy(dirfileres,optionfilefiname);
    strcpy(optfileres,"vpl");
    /* 1eme*/
- Line 8642  set ter svg size 640, 480\nunset log y\n
+ Line 9115  set ter svg size 640, 480\nunset log y\n
              fprintf(ficgp," u %d:(",ioffset);
              kl=0;
              strcpy(gplotcondition,"(");
-             for (k=1; k<=cptcoveff; k++){    /* For each covariate writing the chain of conditions */
+             /* for (k=1; k<=cptcoveff; k++){    /\* For each covariate writing the chain of conditions *\/ */
                /* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate value corresponding to combination k1 and covariate k *\/ */
-               lv=codtabm(k1,TnsdVar[Tvaraff[k]]);
+             for (k=1; k<=cptcovs; k++){    /* For each covariate k get corresponding value lv for combination k1 */
+               /* lv=codtabm(k1,TnsdVar[Tvaraff[k]]); */
+               lv=Tvresult[nres][k];
+               vlv=TinvDoQresult[nres][Tvresult[nres][k]];
                /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
                /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
                /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
                /* vlv= nbcode[Tvaraff[k]][lv]; /\* Value of the modality of Tvaraff[k] *\/ */
-               vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
+               /* vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])]; */
                kl++;
-               sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]);
+               /* sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); */
+               sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,lv, kl+1, vlv );
                kl++;
-               if(k <cptcoveff && cptcoveff>1)
+               if(k <cptcovs && cptcovs>1)
                  sprintf(gplotcondition+strlen(gplotcondition)," && ");
              }
              strcpy(gplotcondition+strlen(gplotcondition),")");
- Line 8737  set ter svg size 640, 480\nunset log y\n
+ Line 9214  set ter svg size 640, 480\nunset log y\n
            }else{
              fprintf(ficgp,",\\\n '' ");
            }
-           if(cptcoveff ==0){ /* No covariate */
+           /* if(cptcoveff ==0){ /\* No covariate *\/ */
+           if(cptcovs ==0){ /* No covariate */
              ioffset=2; /* Age is in 2 */
              /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
              /*#   1       2   3   4   5  6    7  8   9   10  11  12  13  14  15  16  17  18 */
- Line 8849  set ter svg size 640, 480\nunset log y\n
+ Line 9327  set ter svg size 640, 480\nunset log y\n
      fprintf(ficgp,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year\n");
      fprintf(ficgp,"#model=1+age+%s \n",model);
      fprintf(ficgp,"# Type of graphic ng=%d\n",ng);
-     fprintf(ficgp,"#   k1=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */
+     /* fprintf(ficgp,"#   k1=1 to 2^%d=%d\n",cptcoveff,m);/\* to be checked *\/ */
+     fprintf(ficgp,"#   k1=1 to 2^%d=%d\n",cptcovs,m);/* to be checked */
      /* for(k1=1; k1 <=m; k1++)  /\* For each combination of covariate *\/ */
      for(nres=1; nres <= nresult; nres++){ /* For each resultline */
       /* k1=nres; */
- Line 8974  set ter svg size 640, 480\nunset log y\n
+ Line 9453  set ter svg size 640, 480\nunset log y\n
                    } /* end Tprod */
                  }
                  break;
+               case 3:
+                 if(cptcovdageprod >0){
+                   /* if(j==Tprod[ijp]) { */ /* not necessary */
+                     /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
+                     if(ijp <=cptcovprod) { /* Product Vn*Vm and age*VN*Vm*/
+                       if(DummyV[Tvardk[ijp][1]]==0){/* Vn is dummy */
+                         if(DummyV[Tvardk[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*x",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*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvardk[ijp][1]],Tqinvresult[nres][Tvardk[ijp][2]]);
+                         }
+                       }else{ /* age* Vn*Vm Vn is quanti HERE */
+                         if(DummyV[Tvard[ijp][2]]==0){
+                           fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvardk[ijp][2]],Tqinvresult[nres][Tvardk[ijp][1]]);
+                         }else{ /* Both quanti */
+                           fprintf(ficgp,"+p%d*%f*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvardk[ijp][1]],Tqinvresult[nres][Tvardk[ijp][2]]);
+                         }
+                       }
+                       ijp++;
+                     }
+                     /* } */ /* 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 *\/ */
- Line 9060  set ter svg size 640, 480\nunset log y\n
+ Line 9564  set ter svg size 640, 480\nunset log y\n
                      } /* end Tprod */
                    } /* end if */
                    break;
+                 case 3:
+                   if(cptcovdageprod >0){
+                     /* if(j==Tprod[ijp]) { /\* *\/  */
+                       /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
+                       if(ijp <=cptcovprod) { /* Product */
+                         if(DummyV[Tvardk[ijp][1]]==0){/* Vn is dummy */
+                           if(DummyV[Tvardk[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*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvardk[ijp][1]],Tinvresult[nres][Tvardk[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*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvardk[ijp][1]],Tqinvresult[nres][Tvardk[ijp][2]]);
+                             /* fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); */
+                           }
+                         }else{ /* Vn*Vm Vn is quanti */
+                           if(DummyV[Tvardk[ijp][2]]==0){
+                             fprintf(ficgp,"+p%d*%d*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvardk[ijp][2]],Tqinvresult[nres][Tvardk[ijp][1]]);
+                             /* fprintf(ficgp,"+p%d*%d*%f*x",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][Tvardk[ijp][1]],Tqinvresult[nres][Tvardk[ijp][2]]);
+                             /* fprintf(ficgp,"+p%d*%f*%f*x",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 *\/ */
- Line 9413  void prevforecast(char fileres[], double
+ Line 9946  void prevforecast(char fileres[], double
    /* date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); */
    /* date2dmy(dateprojd,&jprojd, &mprojd, &anprojd); */
    /* date2dmy(dateprojf,&jprojf, &mprojf, &anprojf); */
-   i1=pow(2,cptcoveff);
+   /* i1=pow(2,cptcoveff); */
-   if (cptcovn < 1){i1=1;}
+   /* if (cptcovn < 1){i1=1;} */
    fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2);
    fprintf(ficresf,"#****** Routine prevforecast **\n");
  /*            if (h==(int)(YEARM*yearp)){ */
-   for(nres=1; nres <= nresult; nres++) /* For each resultline */
+   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-     for(k=1; k<=i1;k++){ /* We want to find the combination k corresponding to the values of the dummies given in this resut line (to be cleaned one day) */
+     k=TKresult[nres];
-     if(i1 != 1 && TKresult[nres]!= k)
+     if(TKresult[nres]==0) k=1; /* To be checked for noresult */
-       continue;
+     /*  for(k=1; k<=i1;k++){ /\* We want to find the combination k corresponding to the values of the dummies given in this resut line (to be cleaned one day) *\/ */
-     if(invalidvarcomb[k]){
+     /* if(i1 != 1 && TKresult[nres]!= k) */
-       printf("\nCombination (%d) projection ignored because no cases \n",k);
+     /*   continue; */
-       continue;
+     /* if(invalidvarcomb[k]){ */
-     }
+     /*   printf("\nCombination (%d) projection ignored because no cases \n",k);  */
+     /*   continue; */
+     /* } */
      fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#");
-     for(j=1;j<=cptcoveff;j++) {
+     for(j=1;j<=cptcovs;j++){
-       /* fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); */
+       /* for(j=1;j<=cptcoveff;j++) { */
-       fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+     /*   /\* fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); *\/ */
-     }
+     /*   fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[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]);
+     /* for (k4=1; k4<= nsq; k4++){ /\* For each selected (single) quantitative value *\/ */
+     /*   fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */
+     /* } */
+       fprintf(ficresf," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
      }
      fprintf(ficresf," yearproj age");
      for(j=1; j<=nlstate+ndeath;j++){
        for(i=1; i<=nlstate;i++)
- Line 9461  void prevforecast(char fileres[], double
+ Line 10000  void prevforecast(char fileres[], double
            }
          }
          fprintf(ficresf,"\n");
-         for(j=1;j<=cptcoveff;j++)
+         /* for(j=1;j<=cptcoveff;j++)  */
+         for(j=1;j<=cptcovs;j++)
+           fprintf(ficresf,"%d %lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
            /* fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); /\* Tvaraff not correct *\/ */
-           fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); /* TnsdVar[Tvaraff]  correct */
+           /* fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); /\* TnsdVar[Tvaraff]  correct *\/ */
          fprintf(ficresf,"%.f %.f ",anprojd+yearp,agec+h*hstepm/YEARM*stepm);
          for(j=1; j<=nlstate+ndeath;j++) {
- Line 9555  void prevforecast(char fileres[], double
+ Line 10096  void prevforecast(char fileres[], double
    /* if(jintmean==0) jintmean=1; */
    /* if(mintmean==0) jintmean=1; */
-   i1=pow(2,cptcoveff);
+   /* i1=pow(2,cptcoveff); */
-   if (cptcovn < 1){i1=1;}
+   /* if (cptcovn < 1){i1=1;} */
    fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2);
    printf("# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2);
    fprintf(ficresfb,"#****** Routine prevbackforecast **\n");
-   for(nres=1; nres <= nresult; nres++) /* For each resultline */
+   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-   for(k=1; k<=i1;k++){
+     k=TKresult[nres];
-     if(i1 != 1 && TKresult[nres]!= k)
+     if(TKresult[nres]==0) k=1; /* To be checked for noresult */
-       continue;
+   /* for(k=1; k<=i1;k++){ */
-     if(invalidvarcomb[k]){
+   /*   if(i1 != 1 && TKresult[nres]!= k) */
-       printf("\nCombination (%d) projection ignored because no cases \n",k);
+   /*     continue; */
-       continue;
+   /*   if(invalidvarcomb[k]){ */
-     }
+   /*     printf("\nCombination (%d) projection ignored because no cases \n",k);  */
+   /*     continue; */
+   /*   } */
      fprintf(ficresfb,"\n#****** hbijx=probability over h years, hb.jx is weighted by observed prev \n#");
-     for(j=1;j<=cptcoveff;j++) {
+     for(j=1;j<=cptcovs;j++){
-       fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+     /* for(j=1;j<=cptcoveff;j++) { */
-     }
+     /*   fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[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]);
+       fprintf(ficresfb," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
      }
+    /*  fprintf(ficrespij,"******\n"); */
+    /* for (k4=1; k4<= nsq; k4++){ /\* For each selected (single) quantitative value *\/ */
+    /*    fprintf(ficresfb," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */
+    /*  } */
      fprintf(ficresfb," yearbproj age");
      for(j=1; j<=nlstate+ndeath;j++){
        for(i=1; i<=nlstate;i++)
- Line 9608  void prevforecast(char fileres[], double
+ Line 10155  void prevforecast(char fileres[], double
            }
          }
          fprintf(ficresfb,"\n");
-         for(j=1;j<=cptcoveff;j++)
+         /* for(j=1;j<=cptcoveff;j++) */
-           fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+         for(j=1;j<=cptcovs;j++)
+           fprintf(ficresfb,"%d %lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
+           /* fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
          fprintf(ficresfb,"%.f %.f ",anbackd+yearp,agec-h*hstepm/YEARM*stepm);
          for(i=1; i<=nlstate+ndeath;i++) {
            ppij=0.;ppi=0.;
- Line 10210  int readdata(char datafile[], int firsto
+ Line 10759  int readdata(char datafile[], int firsto
    char stra[MAXLINE], strb[MAXLINE];
    char *stratrunc;
-   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++){
-     DummyV[v]=0;
-     FixedV[v]=0;
-   }
-   for(v=ncovcol+1; v <=ncovcol+nqv;v++){
-     DummyV[v]=1;
-     FixedV[v]=0;
-   }
-   for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){
-     DummyV[v]=0;
-     FixedV[v]=1;
-   }
-   for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){
-     DummyV[v]=1;
-     FixedV[v]=1;
-   }
-   for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){
-     printf("Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
-     fprintf(ficlog,"Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
-   }
    ncovcolt=ncovcol+nqv+ntv+nqtv; /* total of covariates in the data, not in the model equation */
- Line 10317  int readdata(char datafile[], int firsto
+ Line 10841  int readdata(char datafile[], int firsto
          if(strb[0]=='.') { /* Missing value */
            lval=-1;
            cotqvar[j][iv][i]=-1; /* 0.0/0.0 */
-           cotvar[j][ntv+iv][i]=-1; /* For performance reasons */
+           cotvar[j][ncovcol+nqv+ntv+iv][i]=-1; /* For performance reasons */
            if(isalpha(strb[1])) { /* .m or .d Really Missing value */
              printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);
              fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value.  Exiting.\n", strb, linei,i,line,iv, nqtv, j);fflush(ficlog);
- Line 10337  int readdata(char datafile[], int firsto
+ Line 10861  int readdata(char datafile[], int firsto
              return 1;
            }
            cotqvar[j][iv][i]=dval;
-           cotvar[j][ntv+iv][i]=dval;
+           cotvar[j][ncovcol+nqv+ntv+iv][i]=dval; /* because cotvar starts now at first ntv */
          }
          strcpy(line,stra);
        }/* end loop ntqv */
- Line 10377  int readdata(char datafile[], int firsto
+ Line 10901  int readdata(char datafile[], int firsto
   Exiting.\n",lval,linei, i,line,iv,j);fflush(ficlog);
            return 1;
          }
-         cotvar[j][iv][i]=(double)(lval);
+         cotvar[j][ncovcol+nqv+iv][i]=(double)(lval);
          strcpy(line,stra);
        }/* end loop ntv */
- Line 10389  int readdata(char datafile[], int firsto
+ Line 10913  int readdata(char datafile[], int firsto
          errno=0;
          lval=strtol(strb,&endptr,10);
          /*      if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
-         if( strb[0]=='\0' || (*endptr != '\0')){
+         if( strb[0]=='\0' || (*endptr != '\0' )){
-           printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);
+           printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);
-           fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
+           fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
+           return 1;
+         }else if( lval==0 || lval > nlstate+ndeath){
+           printf("Error in data around '%s' at line number %d for individual %d, '%s'\n Should be a state at wave %d. A state should be 1 to %d and not %ld.\n Fix your data file '%s'!  Exiting.\n", strb, linei,i,line,j,nlstate+ndeath, lval, datafile);fflush(stdout);
+           fprintf(ficlog,"Error in data around '%s' at line number %d for individual %d, '%s'\n Should be a state at wave %d. A state should be 1 to %d and not %ld.\n Fix your data file '%s'!  Exiting.\n", strb, linei,i,line,j,nlstate+ndeath, lval, datafile); fflush(ficlog);
            return 1;
          }
        }
- Line 10590  int decoderesult( char resultline[], int
+ Line 11118  int decoderesult( char resultline[], int
    printf("decoderesult:%s\n",resultline);
    strcpy(resultsav,resultline);
-   printf("Decoderesult resultsav=\"%s\" resultline=\"%s\"\n", resultsav, resultline);
+   /* printf("Decoderesult resultsav=\"%s\" resultline=\"%s\"\n", resultsav, resultline); */
    if (strlen(resultsav) >1){
      j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' in this resultline */
    }
- Line 10650  int decoderesult( char resultline[], int
+ Line 11178  int decoderesult( char resultline[], int
          if(Tvar[k1]==Tvarsel[k2]) {/* Tvar is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5   */
            modelresult[nres][k2]=k1;/* we found a Vn=1 corrresponding to Vn*age in the model modelresult[2]=1 modelresult[1]=2  modelresult[3]=3  modelresult[6]=4 modelresult[9]=5 */
            resultmodel[nres][k1]=k2; /* Added here */
-           printf("Decoderesult first modelresult[k2=%d]=%d (k1) V%d*AGE\n",k2,k1,Tvar[k1]);
+           /* printf("Decoderesult first modelresult[k2=%d]=%d (k1) V%d*AGE\n",k2,k1,Tvar[k1]); */
            match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */
            break;
          }
- Line 10660  int decoderesult( char resultline[], int
+ Line 11188  int decoderesult( char resultline[], int
          fprintf(ficlog,"Error in result line (Product with age): V%d is missing in result: %s according to model=1+age+%s (Tvarsel[k2=%d]=%d)\n",Tvar[k1], resultline, model, k2, Tvarsel[k2]);
        return 1;
        }
-     }else if(Typevar[k1]==2){ /* Product No age We want to get the position in the resultline of the product in the model line*/
+     }else if(Typevar[k1]==2 || Typevar[k1]==3){ /* Product with or without age. We want to get the position in the resultline of the product in the model line*/
        /* resultmodel[nres][of such a Vn * Vm product k1] is not unique, so can't exist, we feed Tvard[k1][1] and [2] */
        match=0;
-       printf("Decoderesult very first Product Tvardk[k1=%d][1]=%d Tvardk[k1=%d][2]=%d V%d * V%d\n",k1,Tvardk[k1][1],k1,Tvardk[k1][2],Tvardk[k1][1],Tvardk[k1][2]);
+       /* printf("Decoderesult very first Product Tvardk[k1=%d][1]=%d Tvardk[k1=%d][2]=%d V%d * V%d\n",k1,Tvardk[k1][1],k1,Tvardk[k1][2],Tvardk[k1][1],Tvardk[k1][2]); */
        for(k2=1; k2 <=j;k2++){/* Loop on resultline. In result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
          if(Tvardk[k1][1]==Tvarsel[k2]) {/* Tvardk is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5   */
            /* modelresult[k2]=k1; */
-           printf("Decoderesult first Product modelresult[k2=%d]=%d (k1) V%d * \n",k2,k1,Tvarsel[k2]);
+           /* printf("Decoderesult first Product modelresult[k2=%d]=%d (k1) V%d * \n",k2,k1,Tvarsel[k2]); */
            match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */
          }
        }
        if(match == 0){
-         printf("Error in result line (Product without age first variable): V%d is missing in result: %s according to model=1+age+%s\n",Tvardk[k1][1], resultline, model);
+         printf("Error in result line (Product without age first variable or double product with age): V%d is missing in result: %s according to model=1+age+%s\n",Tvardk[k1][1], resultline, model);
-         fprintf(ficlog,"Error in result line (Product without age first variable): V%d is missing in result: %s according to model=1+age+%s\n",Tvardk[k1][1], resultline, model);
+         fprintf(ficlog,"Error in result line (Product without age first variable or double product with age): V%d is missing in result: %s according to model=1+age+%s\n",Tvardk[k1][1], resultline, model);
          return 1;
        }
        match=0;
        for(k2=1; k2 <=j;k2++){/* Loop on resultline. In result line V4=1 V5=24.1 V3=1  V2=8 V1=0 */
          if(Tvardk[k1][2]==Tvarsel[k2]) {/* Tvardk is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5   */
            /* modelresult[k2]=k1;*/
-           printf("Decoderesult second Product modelresult[k2=%d]=%d (k1) * V%d \n ",k2,k1,Tvarsel[k2]);
+           /* printf("Decoderesult second Product modelresult[k2=%d]=%d (k1) * V%d \n ",k2,k1,Tvarsel[k2]); */
            match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */
            break;
          }
        }
        if(match == 0){
-         printf("Error in result line (Product without age second variable): V%d is missing in result: %s according to model=1+age+%s\n",Tvardk[k1][2], resultline, model);
+         printf("Error in result line (Product without age second variable or double product with age): V%d is missing in result: %s according to model=1+age+%s\n",Tvardk[k1][2], resultline, model);
-         fprintf(ficlog,"Error in result line (Product without age second variable): V%d is missing in result : %s according to model=1+age+%s\n",Tvardk[k1][2], resultline, model);
+         fprintf(ficlog,"Error in result line (Product without age second variable or double product with age): V%d is missing in result : %s according to model=1+age+%s\n",Tvardk[k1][2], resultline, model);
          return 1;
        }
      }/* End of testing */
- Line 10699  int decoderesult( char resultline[], int
+ Line 11227  int decoderesult( char resultline[], int
      match=0;
      for(k1=1; k1<= cptcovt ;k1++){ /* loop on model: model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
        if(Typevar[k1]==0){ /* Single only */
-         if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=4 == Tvarsel[1]=4   */
+         if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=4 == Tvarsel[1]=4  What if a product?  */
            resultmodel[nres][k1]=k2;  /* k1th position in the model equation corresponds to k2th position in the result line. resultmodel[2]=1 resultmodel[1]=2  resultmodel[3]=3  resultmodel[6]=4 resultmodel[9]=5 */
            modelresult[nres][k2]=k1; /* k1th position in the model equation corresponds to k2th position in the result line. modelresult[1]=2 modelresult[2]=1  modelresult[3]=3  remodelresult[4]=6 modelresult[5]=9 */
            ++match;
- Line 10740  int decoderesult( char resultline[], int
+ Line 11268  int decoderesult( char resultline[], int
      /* 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, k1 is sorting according to MODEL, but k3 to resultline */
+     if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Dummy and Single, fixed or timevarying, k1 is sorting according to MODEL, but k3 to resultline */
         /* k4+1= (not always if quant in model) 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 corresponds to the k1 position in the model line (doesn't work with products)*/
        /* Value in the (current nres) resultline of the variable at the k1th position in the model equation resultmodel[nres][k1]= k3 */
- Line 10764  int decoderesult( char resultline[], int
+ Line 11292  int decoderesult( char resultline[], int
        Tvresult[nres][k3]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */
        Tinvresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */
        precov[nres][k1]=Tvalsel[k3]; /* Value from resultline of the variable at the k1 position in the model */
-       printf("Decoderesult Dummy k=%d, k1=%d precov[nres=%d][k1=%d]=%.f V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k1, nres, k1,precov[nres][k1], k2, k3, (int)Tvalsel[k3], k4);
+       /* printf("Decoderesult Dummy k=%d, k1=%d precov[nres=%d][k1=%d]=%.f V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k1, nres, k1,precov[nres][k1], k2, k3, (int)Tvalsel[k3], k4); */
        k4++;;
      }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                                 */
- Line 10782  int decoderesult( char resultline[], int
+ Line 11310  int decoderesult( char resultline[], int
        Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */
        TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */
        precov[nres][k1]=Tvalsel[k3q];
-       printf("Decoderesult Quantitative nres=%d,precov[nres=%d][k1=%d]=%.f V(k2q=V%d)= Tvalsel[%d]=%d, Tvarsel[%d]=%f\n",nres, nres, k1,precov[nres][k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]);
+       /* printf("Decoderesult Quantitative nres=%d,precov[nres=%d][k1=%d]=%.f V(k2q=V%d)= Tvalsel[%d]=%d, Tvarsel[%d]=%f\n",nres, nres, k1,precov[nres][k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]); */
        k4q++;;
-     }else if( Dummy[k1]==2 ){ /* For dummy with age product */
+     }else if( Dummy[k1]==2 ){ /* For dummy with age product "V2+V3+V4+V6+V7+V6*V2+V7*V2+V6*V3+V7*V3+V6*V4+V7*V4+age*V2+age*V3+age*V4+age*V6+age*V7+age*V6*V2+age*V6*V3+age*V7*V3+age*V6*V4+age*V7*V4\r"*/
-       /* Tvar[k1]; */ /* Age variable */
+       /* Tvar[k1]; */ /* Age variable */ /* 17 age*V6*V2 ?*/
        /* Wrong we want the value of variable name Tvar[k1] */
+       if(Typevar[k1]==2 || Typevar[k1]==3 ){ /* For product quant or dummy (with or without age) */
-       k3= resultmodel[nres][k1]; /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/
+         precov[nres][k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]];
-       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)*/
+       /* printf("Decoderesult Quantitative or Dummy (not with age) nres=%d k1=%d precov[nres=%d][k1=%d]=%.f V%d(=%.f) * V%d(=%.f) \n",nres, k1, nres, k1,precov[nres][k1], Tvardk[k1][1], TinvDoQresult[nres][Tvardk[k1][1]], Tvardk[k1][2], TinvDoQresult[nres][Tvardk[k1][2]]); */
-       TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* TinvDoQresult[nres][4]=1 */
+       }else{
-       precov[nres][k1]=Tvalsel[k3];
+         k3= resultmodel[nres][k1]; /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/
-       printf("Decoderesult Dummy with age k=%d, k1=%d precov[nres=%d][k1=%d]=%.f Tvar[%d]=V%d k2=Tvarsel[%d]=%d Tvalsel[%d]=%d\n",k, k1, nres, k1,precov[nres][k1], k1, Tvar[k1], k3,(int)Tvarsel[k3], k3, (int)Tvalsel[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]; /* TinvDoQresult[nres][4]=1 */
+         precov[nres][k1]=Tvalsel[k3];
+       }
+       /* printf("Decoderesult Dummy with age k=%d, k1=%d precov[nres=%d][k1=%d]=%.f Tvar[%d]=V%d k2=Tvarsel[%d]=%d Tvalsel[%d]=%d\n",k, k1, nres, k1,precov[nres][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][k1]; /* resultmodel[1(V5)] = 25.1=k3q */
+       if(Typevar[k1]==2 || Typevar[k1]==3 ){ /* For product quant or dummy (with or without age) */
-       k2q=(int)Tvarsel[k3q]; /*  Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */
+         precov[nres][k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]];
-       TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* TinvDoQresult[nres][5]=25.1 */
+       /* printf("Decoderesult Quantitative or Dummy (not with age) nres=%d k1=%d precov[nres=%d][k1=%d]=%.f V%d(=%.f) * V%d(=%.f) \n",nres, k1, nres, k1,precov[nres][k1], Tvardk[k1][1], TinvDoQresult[nres][Tvardk[k1][1]], Tvardk[k1][2], TinvDoQresult[nres][Tvardk[k1][2]]); */
-       precov[nres][k1]=Tvalsel[k3q];
+       }else{
-       printf("Decoderesult Quantitative with age nres=%d, k1=%d, precov[nres=%d][k1=%d]=%f Tvar[%d]=V%d V(k2q=%d)= Tvarsel[%d]=%d, Tvalsel[%d]=%f\n",nres, k1, nres, k1,precov[nres][k1], k1,  Tvar[k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]);
+         k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 25.1=k3q */
-     }else if(Typevar[k1]==2 ){ /* For product quant or dummy (not with age) */
+         k2q=(int)Tvarsel[k3q]; /*  Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */
+         TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* TinvDoQresult[nres][5]=25.1 */
+         precov[nres][k1]=Tvalsel[k3q];
+       }
+       /* printf("Decoderesult Quantitative with age nres=%d, k1=%d, precov[nres=%d][k1=%d]=%f Tvar[%d]=V%d V(k2q=%d)= Tvarsel[%d]=%d, Tvalsel[%d]=%f\n",nres, k1, nres, k1,precov[nres][k1], k1,  Tvar[k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]); */
+     }else if(Typevar[k1]==2 || Typevar[k1]==3 ){ /* For product quant or dummy (with or without age) */
        precov[nres][k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]];
-       printf("Decoderesult Quantitative or Dummy (not with age) nres=%d k1=%d precov[nres=%d][k1=%d]=%.f V%d(=%.f) * V%d(=%.f) \n",nres, k1, nres, k1,precov[nres][k1], Tvardk[k1][1], TinvDoQresult[nres][Tvardk[k1][1]], Tvardk[k1][2], TinvDoQresult[nres][Tvardk[k1][2]]);
+       /* printf("Decoderesult Quantitative or Dummy (not with age) nres=%d k1=%d precov[nres=%d][k1=%d]=%.f V%d(=%.f) * V%d(=%.f) \n",nres, k1, nres, k1,precov[nres][k1], Tvardk[k1][1], TinvDoQresult[nres][Tvardk[k1][1]], Tvardk[k1][2], TinvDoQresult[nres][Tvardk[k1][2]]); */
      }else{
        printf("Error Decoderesult probably a product  Dummy[%d]==%d && Typevar[%d]==%d\n", k1, Dummy[k1], k1, Typevar[k1]);
        fprintf(ficlog,"Error Decoderesult probably a product  Dummy[%d]==%d && Typevar[%d]==%d\n", k1, Dummy[k1], k1, Typevar[k1]);
- Line 10831  int decodemodel( char model[], int lasto
+ Line 11368  int decodemodel( char model[], int lasto
  /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */
  {
    int i, j, k, ks, v;
-   int  j1, k1, k2, k3, k4;
+   int n,m;
-   char modelsav[80];
+   int  j1, k1, k11, k12, k2, k3, k4;
-   char stra[80], strb[80], strc[80], strd[80],stre[80];
+   char modelsav[300];
+   char stra[300], strb[300], strc[300], strd[300],stre[300],strf[300];
    char *strpt;
+   int  **existcomb;
+   existcomb=imatrix(1,NCOVMAX,1,NCOVMAX);
+   for(i=1;i<=NCOVMAX;i++)
+     for(j=1;j<=NCOVMAX;j++)
+       existcomb[i][j]=0;
    /*removespace(model);*/
    if (strlen(model) >1){ /* If there is at least 1 covariate */
-     j=0, j1=0, k1=0, k2=-1, ks=0, cptcovn=0;
+     j=0, j1=0, k1=0, k12=0, k2=-1, ks=0, cptcovn=0;
      if (strstr(model,"AGE") !=0){
        printf("Error. AGE must be in lower case 'age' model=1+age+%s. ",model);
        fprintf(ficlog,"Error. AGE must be in lower case model=1+age+%s. ",model);fflush(ficlog);
- Line 10870  int decodemodel( char model[], int lasto
+ Line 11414  int decodemodel( char model[], int lasto
          substrchaine(modelsav, model, "age*age");
      }else
        nagesqr=0;
-     if (strlen(modelsav) >1){
+     if (strlen(modelsav) >1){ /* V2 +V3 +V4 +V6 +V7 +V6*V2 +V7*V2 +V6*V3 +V7*V3 +V6*V4 +V7*V4 +age*V2 +age*V3 +age*V4 +age*V6 +age*V7 +age*V6*V2 +V7*V2 +age*V6*V3 +age*V7*V3 +age*V6*V4 +age*V7*V4 */
        j=nbocc(modelsav,'+'); /**< j=Number of '+' */
        j1=nbocc(modelsav,'*'); /**< j1=Number of '*' */
-       cptcovs=j+1-j1; /**<  Number of simple covariates V1+V1*age+V3 +V3*V4+age*age=> V1 + V3 =5-3=2  */
+       cptcovs=0; /**<  Number of simple covariates V1 +V1*age +V3 +V3*V4 +age*age => V1 + V3 =4+1-3=2  Wrong */
        cptcovt= j+1; /* Number of total covariates in the model, not including
                       * cst, age and age*age
                       * V1+V1*age+ V3 + V3*V4+age*age=> 3+1=4*/
        /* including age products which are counted in cptcovage.
         * but the covariates which are products must be treated
         * separately: ncovn=4- 2=2 (V1+V3). */
-       cptcovprod=j1; /**< Number of products  V1*V2 +v3*age = 2 */
+       cptcovprod=0; /**< Number of products  V1*V2 +v3*age = 2 */
+       cptcovdageprod=0; /* Number of doouble products with age age*Vn*VM or Vn*age*Vm or Vn*Vm*age */
        cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1  */
+       cptcovprodage=0;
+       /* cptcovprodage=nboccstr(modelsav,"age");*/
        /*   Design
         *  V1   V2   V3   V4  V5  V6  V7  V8  V9 Weight
- Line 10890  int decodemodel( char model[], int lasto
+ Line 11436  int decodemodel( char model[], int lasto
         * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8
         *   k=  1    2      3       4     5       6      7        8
         *  cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8
-        *  covar[k,i], value of kth covariate if not including age for individual i:
+        *  covar[k,i], are for fixed covariates, value of kth covariate if not including age for individual i:
         *       covar[1][i]= (V1), covar[4][i]=(V4), covar[8][i]=(V8)
         *  Tvar[k] # of the kth covariate:  Tvar[1]=2  Tvar[2]=1 Tvar[4]=3 Tvar[8]=8
         *       if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and
         *  Tage[++cptcovage]=k
-        *       if products, new covar are created after ncovcol with k1
+        *       if products, new covar are created after ncovcol + nqv (quanti fixed) with k1
         *  Tvar[k]=ncovcol+k1; # of the kth covariate product:  Tvar[5]=ncovcol+1=10  Tvar[6]=ncovcol+1=11
         *  Tprod[k1]=k; Tprod[1]=5 Tprod[2]= 6; gives the position of the k1th product
         *  Tvard[k1][1]=m Tvard[k1][2]=m; Tvard[1][1]=5 (V5) Tvard[1][2]=6 Tvard[2][1]=7 (V7) Tvard[2][2]=8
         *  Tvar[cptcovn+k2]=Tvard[k1][1];Tvar[cptcovn+k2+1]=Tvard[k1][2];
         *  Tvar[8+1]=5;Tvar[8+2]=6;Tvar[8+3]=7;Tvar[8+4]=8 inverted
         *  V1   V2   V3   V4  V5  V6  V7  V8  V9  V10  V11
-        *  <          ncovcol=8                >
+        *  <          ncovcol=8  8 fixed covariate. Additional starts at 9 (V5*V6) and 10(V7*V8)              >
         *       Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8    d1   d1   d2  d2
         *          k=  1    2      3       4     5       6      7        8    9   10   11  12
-        *     Tvar[k]= 2    1      3       3    10      11      8        8    5    6    7   8
+        *     Tvard[k]= 2    1      3       3    10      11      8        8    5    6    7   8
-        * p Tvar[1]@12={2,   1,     3,      3,  11,     10,     8,       8,   7,   8,   5,  6}
+        * p Tvar[1]@12={2,   1,     3,      3,   9,     10,     8,       8}
         * p Tprod[1]@2={                         6, 5}
         *p Tvard[1][1]@4= {7, 8, 5, 6}
         * covar[k][i]= V2   V1      ?      V3    V5*V6?   V7*V8?  ?       V8
- Line 10936  int decodemodel( char model[], int lasto
+ Line 11482  int decodemodel( char model[], int lasto
          Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0;
        }
        cptcovage=0;
+       /* First loop in order to calculate */
+       /* for age*VN*Vm
+        * Provides, Typevar[k], Tage[cptcovage], existcomb[n][m], FixedV[ncovcolt+k12]
+        * Tprod[k1]=k  Tposprod[k]=k1;    Tvard[k1][1] =m;
+       */
+       /* Needs  FixedV[Tvardk[k][1]] */
+       /* For others:
+        * Sets   Typevar[k];
+        * Tvar[k]=ncovcol+nqv+ntv+nqtv+k11;
+        *        Tposprod[k]=k11;
+        *        Tprod[k11]=k;
+        *        Tvardk[k][1] =m;
+        * Needs FixedV[Tvardk[k][1]] == 0
+       */
        for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model line */
          cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' cutl from left to right
                                           modelsav==V2+V1+V5*age+V4+V3*age strb=V3*age stra=V2+V1V5*age+V4 */    /* <model> "V5+V4+V3+V4*V3+V5*age+V1*age+V1" strb="V5" stra="V4+V3+V4*V3+V5*age+V1*age+V1" */
- Line 10943  int decodemodel( char model[], int lasto
+ Line 11505  int decodemodel( char model[], int lasto
            strcpy(strb,modelsav); /* and analyzes it */
          /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
          /*scanf("%d",i);*/
-         if (strchr(strb,'*')) {  /**< Model includes a product V2+V1+V5*age+ V4+V3*age strb=V3*age */
+         if (strchr(strb,'*')) {  /**< Model includes a product V2+V1+V5*age+ V4+V3*age strb=V3*age OR double product with age strb=age*V6*V2 or V6*V2*age or V6*age*V2 */
-           cutl(strc,strd,strb,'*'); /**< k=1 strd*strc  Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
+           cutl(strc,strd,strb,'*'); /**< k=1 strd*strc  Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 OR strb=age*V6*V2 strc=V6*V2 strd=age OR c=V2*age OR c=age*V2  */
-           if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
+           if(strchr(strc,'*')) { /**< Model with age and DOUBLE product: allowed since 0.99r44, strc=V6*V2 or V2*age or age*V2, strd=age or V6 or V6   */
-             /* covar is not filled and then is empty */
+             Typevar[k]=3;  /* 3 for age and double product age*Vn*Vm varying of fixed */
-             cptcovprod--;
+             if(strstr(strc,"age")!=0) { /* It means that strc=V2*age or age*V2 and thus that strd=Vn */
-             cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
+               cutl(stre,strf,strc,'*') ; /* strf=age or Vm, stre=Vm or age. If strc=V6*V2 then strf=V6 and stre=V2 */
-             Tvar[k]=atoi(stre);  /* V2+V1+V5*age+V4+V3*age Tvar[5]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */
+               strcpy(strc,strb); /* save strb(=age*Vn*Vm) into strc */
-             Typevar[k]=1;  /* 1 for age product */
+               /* We want strb=Vn*Vm */
-             cptcovage++; /* Counts the number of covariates which include age as a product */
+               if(strcmp(strf,"age")==0){ /* strf is "age" so that stre=Vm =V2 . */
-             Tage[cptcovage]=k;  /*  V2+V1+V4+V3*age Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
+                 strcpy(strb,strd);
-             /*printf("stre=%s ", stre);*/
+                 strcat(strb,"*");
-           } else if (strcmp(strd,"age")==0) { /* or age*Vn */
+                 strcat(strb,stre);
-             cptcovprod--;
+               }else{  /* strf=Vm  If strf=V6 then stre=V2 */
-             cutl(stre,strb,strc,'V');
+                 strcpy(strb,strf);
-             Tvar[k]=atoi(stre);
+                 strcat(strb,"*");
-             Typevar[k]=1;  /* 1 for age product */
+                 strcat(strb,stre);
-             cptcovage++;
+                 strcpy(strd,strb); /* in order for strd to not be "age"  for next test (will be Vn*Vm */
-             Tage[cptcovage]=k;
+               }
-           } else {  /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2  strb=V3*V2*/
+               /* printf("DEBUG FIXED k=%d, Tage[k]=%d, Tvar[Tage[k]=%d,FixedV[Tvar[Tage[k]]]=%d\n",k,Tage[k],Tvar[Tage[k]],FixedV[Tvar[Tage[k]]]); */
-             /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
+               /* FixedV[Tvar[Tage[k]]]=0; /\* HERY not sure if V7*V4*age Fixed might not exist  yet*\/ */
-             cptcovn++;
+             }else{  /* strc=Vn*Vm (and strd=age) and should be strb=Vn*Vm but want to keep original strb double product  */
-             cptcovprodnoage++;k1++;
+               strcpy(stre,strb); /* save full b in stre */
+               strcpy(strb,strc); /* save short c in new short b for next block strb=Vn*Vm*/
+               strcpy(strf,strc); /* save short c in new short f */
+               cutl(strc,strd,strf,'*'); /* We get strd=Vn and strc=Vm for next block (strb=Vn*Vm)*/
+               /* strcpy(strc,stre);*/ /* save full e in c for future */
+             }
+             cptcovdageprod++; /* double product with age  Which product is it? */
+             /* strcpy(strb,strc);  /\* strb was age*V6*V2 or V6*V2*age or V6*age*V2 IS now V6*V2 or V2*age or age*V2 *\/ */
+             /* cutl(strc,strd,strb,'*'); /\* strd=  V6    or   V2     or    age and  strc=  V2 or    age or    V2 *\/ */
              cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
-             Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* ncovcolt+k1; For model-covariate k tells which data-covariate to use but
+             n=atoi(stre);
-                                                 because this model-covariate is a construction we invent a new column
-                                                 which is after existing variables ncovcol+nqv+ntv+nqtv + k1
-                                                 If already ncovcol=4 and model= V2 + V1 + V1*V4 + age*V3 + V3*V2
-                                                 thus after V4 we invent V5 and V6 because age*V3 will be computed in 4
-                                                 Tvar[3=V1*V4]=4+1=5 Tvar[5=V3*V2]=4 + 2= 6, Tvar[4=age*V3]=3 etc */
-             /* Please remark that the new variables are model dependent */
-             /* If we have 4 variable but the model uses only 3, like in
-              * model= V1 + age*V1 + V2 + V3 + age*V2 + age*V3 + V1*V2 + V1*V3
-              *  k=     1     2       3   4     5        6        7       8
-              * Tvar[k]=1     1       2   3     2        3       (5       6) (and not 4 5 because of V4 missing)
-              * Tage[kk]    [1]= 2           [2]=5      [3]=6                  kk=1 to cptcovage=3
-              * Tvar[Tage[kk]][1]=2          [2]=2      [3]=3
-              */
-             Typevar[k]=2;  /* 2 for product */
              cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
-             Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */
+             m=atoi(strc);
-             Tposprod[k]=k1; /* Tposprod[3]=1, Tposprod[2]=5 */
+             cptcovage++; /* Counts the number of covariates which include age as a product */
-             Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
+             Tage[cptcovage]=k; /* For age*V3*V2 gives the position in model of covariates associated with age Tage[1]=6 HERY too*/
-             Tvardk[k][1] =atoi(strc); /* m 1 for V1*/
+             if(existcomb[n][m] == 0){
-             Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
+               /*  r /home/brouard/Documents/Recherches/REVES/Zachary/Zach-2022/Feinuo_Sun/Feinuo-threeway/femV12V15_3wayintNBe.imach */
-             Tvardk[k][2] =atoi(stre); /* n 4 for V4*/
+               printf("Warning in model combination V%d*V%d should exist in the model before adding V%d*V%d*age !\n",n,m,n,m);
-             k2=k2+2;  /* k2 is initialize to -1, We want to store the n and m in Vn*Vm at the end of Tvar */
+               fprintf(ficlog,"Warning in model combination V%d*V%d should exist in the model before adding V%d*V%d*age !\n",n,m,n,m);
-             /* Tvar[cptcovt+k2]=Tvard[k1][1]; /\* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) *\/ */
+               fflush(ficlog);
-             /* Tvar[cptcovt+k2+1]=Tvard[k1][2];  /\* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) *\/ */
+               k1++;  /* The combination Vn*Vm will be in the model so we create it at k1 */
-             /*ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2, Tvar[3]=5, Tvar[4]=6, cptcovt=5 */
+               k12++;
-             /*                     1  2   3      4     5 | Tvar[5+1)=1, Tvar[7]=2   */
+               existcomb[n][m]=k1;
-             if( FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){ /* If the product is a fixed covariate then we feed the new column with Vn*Vm */
+               existcomb[m][n]=k1;
-               for (i=1; i<=lastobs;i++){/* For fixed product */
+               Tvar[k]=ncovcol+nqv+ntv+nqtv+k1;
-               /* Computes the new covariate which is a product of
+               Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2+ age*V6*V3 Gives the k position of the k1 double product Vn*Vm or age*Vn*Vm*/
-                  covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
+               Tposprod[k]=k1; /* Tposprod[3]=1, Tposprod[2]=5 Gives the k1 double product  Vn*Vm or age*Vn*Vm at the k position */
-               covar[ncovcolt+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
+               Tvard[k1][1] =m; /* m 1 for V1*/
+               Tvardk[k][1] =m; /* m 1 for V1*/
+               Tvard[k1][2] =n; /* n 4 for V4*/
+               Tvardk[k][2] =n; /* n 4 for V4*/
+ /*            Tvar[Tage[cptcovage]]=k1;*/ /* Tvar[6=age*V3*V2]=9 (new fixed covariate) */ /* We don't know about Fixed yet HERE */
+               if( FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){ /* If the product is a fixed covariate then we feed the new column with Vn*Vm */
+                 for (i=1; i<=lastobs;i++){/* For fixed product */
+                   /* Computes the new covariate which is a product of
+                      covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
+                   covar[ncovcolt+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
+                 }
+                 cptcovprodage++; /* Counting the number of fixed covariate with age */
+                 FixedV[ncovcolt+k12]=0; /* We expand Vn*Vm */
+                 k12++;
+                 FixedV[ncovcolt+k12]=0;
+               }else{ /*End of FixedV */
+                 cptcovprodvage++; /* Counting the number of varying covariate with age */
+                 FixedV[ncovcolt+k12]=1; /* We expand Vn*Vm */
+                 k12++;
+                 FixedV[ncovcolt+k12]=1;
                }
-             } /*End of FixedV */
+             }else{  /* k1 Vn*Vm already exists */
-           } /* End age is not in the model */
+               k11=existcomb[n][m];
-         } /* End if model includes a product */
+               Tposprod[k]=k11; /* OK */
-         else { /* not a product */
+               Tvar[k]=Tvar[Tprod[k11]]; /* HERY */
+               Tvardk[k][1]=m;
+               Tvardk[k][2]=n;
+               if( FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){ /* If the product is a fixed covariate then we feed the new column with Vn*Vm */
+                 /*cptcovage++;*/ /* Counts the number of covariates which include age as a product */
+                 cptcovprodage++; /* Counting the number of fixed covariate with age */
+                 /*Tage[cptcovage]=k;*/ /* For age*V3*V2 Tage[1]=V3*V3=9 HERY too*/
+                 Tvar[Tage[cptcovage]]=k1;
+                 FixedV[ncovcolt+k12]=0; /* We expand Vn*Vm */
+                 k12++;
+                 FixedV[ncovcolt+k12]=0;
+               }else{ /* Already exists but time varying (and age) */
+                 /*cptcovage++;*/ /* Counts the number of covariates which include age as a product */
+                 /*Tage[cptcovage]=k;*/ /* For age*V3*V2 Tage[1]=V3*V3=9 HERY too*/
+                 /* Tvar[Tage[cptcovage]]=k1; */
+                 cptcovprodvage++;
+                 FixedV[ncovcolt+k12]=1; /* We expand Vn*Vm */
+                 k12++;
+                 FixedV[ncovcolt+k12]=1;
+               }
+             }
+             /* Tage[cptcovage]=k;  /\*  V2+V1+V4+V3*age Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 *\/ */
+             /* Tvar[k]=k11; /\* HERY *\/ */
+           } else {/* simple product strb=age*Vn so that c=Vn and d=age, or strb=Vn*age so that c=age and d=Vn, or b=Vn*Vm so that c=Vm and d=Vn */
+             cptcovprod++;
+             if (strcmp(strc,"age")==0) { /**< Model includes age: strb= Vn*age c=age d=Vn*/
+               /* covar is not filled and then is empty */
+               cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
+               Tvar[k]=atoi(stre);  /* V2+V1+V5*age+V4+V3*age Tvar[5]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */
+               Typevar[k]=1;  /* 1 for age product */
+               cptcovage++; /* Counts the number of covariates which include age as a product */
+               Tage[cptcovage]=k;  /*  V2+V1+V4+V3*age Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
+               if( FixedV[Tvar[k]] == 0){
+                 cptcovprodage++; /* Counting the number of fixed covariate with age */
+               }else{
+                 cptcovprodvage++; /* Counting the number of fixedvarying covariate with age */
+               }
+               /*printf("stre=%s ", stre);*/
+             } else if (strcmp(strd,"age")==0) { /* strb= age*Vn c=Vn */
+               cutl(stre,strb,strc,'V');
+               Tvar[k]=atoi(stre);
+               Typevar[k]=1;  /* 1 for age product */
+               cptcovage++;
+               Tage[cptcovage]=k;
+               if( FixedV[Tvar[k]] == 0){
+                 cptcovprodage++; /* Counting the number of fixed covariate with age */
+               }else{
+                 cptcovprodvage++; /* Counting the number of fixedvarying covariate with age */
+               }
+             }else{ /*  for product Vn*Vm */
+               Typevar[k]=2;  /* 2 for product Vn*Vm */
+               cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
+               n=atoi(stre);
+               cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
+               m=atoi(strc);
+               k1++;
+               cptcovprodnoage++;
+               if(existcomb[n][m] != 0 || existcomb[m][n] != 0){
+                 printf("Warning in model combination V%d*V%d already exists in the model in position k1=%d!\n",n,m,existcomb[n][m]);
+                 fprintf(ficlog,"Warning in model combination V%d*V%d already exists in the model in position k1=%d!\n",n,m,existcomb[n][m]);
+                 fflush(ficlog);
+                 k11=existcomb[n][m];
+                 Tvar[k]=ncovcol+nqv+ntv+nqtv+k11;
+                 Tposprod[k]=k11;
+                 Tprod[k11]=k;
+                 Tvardk[k][1] =m; /* m 1 for V1*/
+                 /* Tvard[k11][1] =m; /\* n 4 for V4*\/ */
+                 Tvardk[k][2] =n; /* n 4 for V4*/
+                 /* Tvard[k11][2] =n; /\* n 4 for V4*\/ */
+               }else{ /* combination Vn*Vm doesn't exist we create it (no age)*/
+                 existcomb[n][m]=k1;
+                 existcomb[m][n]=k1;
+                 Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* ncovcolt+k1; For model-covariate k tells which data-covariate to use but
+                                                     because this model-covariate is a construction we invent a new column
+                                                     which is after existing variables ncovcol+nqv+ntv+nqtv + k1
+                                                     If already ncovcol=4 and model= V2 + V1 + V1*V4 + age*V3 + V3*V2
+                                                     thus after V4 we invent V5 and V6 because age*V3 will be computed in 4
+                                                     Tvar[3=V1*V4]=4+1=5 Tvar[5=V3*V2]=4 + 2= 6, Tvar[4=age*V3]=3 etc */
+                 /* Please remark that the new variables are model dependent */
+                 /* If we have 4 variable but the model uses only 3, like in
+                  * model= V1 + age*V1 + V2 + V3 + age*V2 + age*V3 + V1*V2 + V1*V3
+                  *  k=     1     2      3   4     5        6        7       8
+                  * Tvar[k]=1     1       2   3     2        3       (5       6) (and not 4 5 because of V4 missing)
+                  * Tage[kk]    [1]= 2           [2]=5      [3]=6                  kk=1 to cptcovage=3
+                  * Tvar[Tage[kk]][1]=2          [2]=2      [3]=3
+                  */
+                 /* We need to feed some variables like TvarVV, but later on next loop because of ncovv (k2) is not correct */
+                 Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 +V6*V2*age  */
+                 Tposprod[k]=k1; /* Tposprod[3]=1, Tposprod[2]=5 */
+                 Tvard[k1][1] =m; /* m 1 for V1*/
+                 Tvardk[k][1] =m; /* m 1 for V1*/
+                 Tvard[k1][2] =n; /* n 4 for V4*/
+                 Tvardk[k][2] =n; /* 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) *\/ */
+                 /*ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2, Tvar[3]=5, Tvar[4]=6, cptcovt=5 */
+                 /*                     1  2   3      4     5 | Tvar[5+1)=1, Tvar[7]=2   */
+                 if( FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){ /* If the product is a fixed covariate then we feed the new column with Vn*Vm */
+                   for (i=1; i<=lastobs;i++){/* For fixed product */
+                     /* Computes the new covariate which is a product of
+                        covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
+                     covar[ncovcolt+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
+                   }
+                   /* TvarVV[k2]=n; */
+                   /* FixedV[ncovcolt+k2]=0; /\* or FixedV[Tvar[k]]=0; FixedV[TvarVV[ncovv]]=0 HERE *\/ */
+                   /* TvarVV[k2+1]=m; */
+                   /* FixedV[ncovcolt+k2]=0; /\* or FixedV[Tvar[k]]=0; FixedV[TvarVV[ncovv]]=0 HERE *\/ */
+                 }else{ /* not FixedV */
+                   /* TvarVV[k2]=n; */
+                   /* FixedV[ncovcolt+k2]=0; /\* or FixedV[Tvar[k]]=0; FixedV[TvarVV[ncovv]]=0 HERE *\/ */
+                   /* TvarVV[k2+1]=m; */
+                   /* FixedV[ncovcolt+k2]=0; /\* or FixedV[Tvar[k]]=0; FixedV[TvarVV[ncovv]]=0 HERE *\/ */
+                 }
+               }  /* End of creation of Vn*Vm if not created by age*Vn*Vm earlier  */
+             } /*  End of product Vn*Vm */
+           } /* End of age*double product or simple product */
+         }else { /* not a product */
            /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
            /*  scanf("%d",i);*/
            cutl(strd,strc,strb,'V');
- Line 11015  int decodemodel( char model[], int lasto
+ Line 11707  int decodemodel( char model[], int lasto
                                  /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
                                    scanf("%d",i);*/
        } /* end of loop + on total covariates */
      } /* end if strlen(modelsave == 0) age*age might exist */
    } /* end if strlen(model == 0) */
+   cptcovs=cptcovt - cptcovdageprod - cptcovprod;/**<  Number of simple covariates V1 +V1*age +V3 +V3*V4 +age*age + age*v4*V3=> V1 + V3 =4+1-3=2  */
    /*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products.
      If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/
- Line 11043  int decodemodel( char model[], int lasto
+ Line 11738  int decodemodel( char model[], int lasto
    /* Tvar[k] is the value n of Vn with n varying for 1 to nvcol, or p  Vp=Vn*Vm for product */
          /* Computing effective variables, ie used by the model, that is from the cptcovt variables */
    printf("Model=1+age+%s\n\
- Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product \n\
+ Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product, 3 for double product with age \n\
  Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
  Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
    fprintf(ficlog,"Model=1+age+%s\n\
- Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product \n\
+ Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product, 3 for double product with age  \n\
  Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
  Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
-   for(k=-1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;}
+   for(k=-1;k<=NCOVMAX; k++){ Fixed[k]=0; Dummy[k]=0;}
-   for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0, ncovvt=0;k<=cptcovt; k++){ /* or cptocvt */
+   for(k=1;k<=NCOVMAX; k++){TvarFind[k]=0; TvarVind[k]=0;}
+   /* Second loop for calculating  Fixed[k], Dummy[k]*/
+   for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0,ncovva=0,ncovvta=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0, ncovvt=0;k<=cptcovt; k++){ /* or cptocvt loop on k from model */
      if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
        Fixed[k]= 0;
        Dummy[k]= 0;
- Line 11067  Dummy[k] 0=dummy (0 1), 1 quantitative (
+ Line 11768  Dummy[k] 0=dummy (0 1), 1 quantitative (
        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 */
      /* }else if( Tvar[k] <=ncovcol &&  Typevar[k]==2){ /\* Product of fixed dummy (<=ncovcol) covariates For a fixed product k is higher than ncovcol *\/ */
-     }else if( Tposprod[k]>0  &&  Typevar[k]==2 && FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){ /* Needs a fixed product Product of fixed dummy (<=ncovcol) covariates For a fixed product k is higher than ncovcol */
-       Fixed[k]= 0;
-       Dummy[k]= 0;
-       ncoveff++;
-       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 */
      }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){/* Remind that product Vn*Vm are added in k Only simple fixed quantitative variable */
        Fixed[k]= 0;
        Dummy[k]= 1;
- Line 11139  Dummy[k] 0=dummy (0 1), 1 quantitative (
+ Line 11829  Dummy[k] 0=dummy (0 1), 1 quantitative (
        TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */
        TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */
        /* Tmodeliqind[k]=nqtveff;/\* Only simple time varying quantitative variable *\/ */
-       printf("Quasi TmodelQind[%d]=%d,Tvar[TmodelQind[%d]]=V%d, ncovcol=%d, nqv=%d, ntv=%d,Tvar[k]- ncovcol-nqv-ntv=%d\n",nqtveff,k,nqtveff,Tvar[k], ncovcol, nqv, ntv, Tvar[k]- ncovcol-nqv-ntv);
+       /* printf("Quasi TmodelQind[%d]=%d,Tvar[TmodelQind[%d]]=V%d, ncovcol=%d, nqv=%d, ntv=%Ad,Tvar[k]- ncovcol-nqv-ntv=%d\n",nqtveff,k,nqtveff,Tvar[k], ncovcol, nqv, ntv, Tvar[k]- ncovcol-nqv-ntv); */
-       printf("Quasi TmodelInvQind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv-ntv);
+       /* printf("Quasi TmodelInvQind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv-ntv); */
      }else if (Typevar[k] == 1) {  /* product with age */
        ncova++;
        TvarA[ncova]=Tvar[k];
        TvarAind[ncova]=k;
+       /** Fixed[k] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */
+       /** Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */
        if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */
          Fixed[k]= 2;
          Dummy[k]= 2;
          modell[k].maintype= ATYPE;
          modell[k].subtype= APFD;
+         ncovta++;
+         TvarAVVA[ncovta]=Tvar[k]; /*  (2)age*V3 */
+         TvarAVVAind[ncovta]=k;
          /* ncoveff++; */
        }else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/
          Fixed[k]= 2;
          Dummy[k]= 3;
          modell[k].maintype= ATYPE;
          modell[k].subtype= APFQ;                /*      Product age * fixed quantitative */
+         ncovta++;
+         TvarAVVA[ncovta]=Tvar[k]; /*   */
+         TvarAVVAind[ncovta]=k;
          /* nqfveff++;  /\* Only simple fixed quantitative variable *\/ */
        }else if( Tvar[k] <=ncovcol+nqv+ntv ){
          Fixed[k]= 3;
          Dummy[k]= 2;
          modell[k].maintype= ATYPE;
          modell[k].subtype= APVD;                /*      Product age * varying dummy */
+         ncovva++;
+         TvarVVA[ncovva]=Tvar[k]; /*  (1)+age*V6 + (2)age*V7 */
+         TvarVVAind[ncovva]=k;
+         ncovta++;
+         TvarAVVA[ncovta]=Tvar[k]; /*   */
+         TvarAVVAind[ncovta]=k;
          /* ntveff++; /\* Only simple time varying dummy variable *\/ */
        }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){
          Fixed[k]= 3;
          Dummy[k]= 3;
          modell[k].maintype= ATYPE;
          modell[k].subtype= APVQ;                /*      Product age * varying quantitative */
+         ncovva++;
+         TvarVVA[ncovva]=Tvar[k]; /*   */
+         TvarVVAind[ncovva]=k;
+         ncovta++;
+         TvarAVVA[ncovta]=Tvar[k]; /*  (1)+age*V6 + (2)age*V7 */
+         TvarAVVAind[ncovta]=k;
          /* nqtveff++;/\* Only simple time varying quantitative variable *\/ */
        }
-     }else if (Typevar[k] == 2) {  /* product Vn * Vm without age, V1+V3+age*V1+age*V3+V1*V3 looking at V1*V3, Typevar={0, 0, 1, 1, 2}, k=5, V1 is fixed, V3 is timevary, V5 is a product  */
+     }else if( Tposprod[k]>0  &&  Typevar[k]==2){  /* Detects if fixed product no age Vm*Vn */
+       printf("MEMORY ERRORR k=%d  Tposprod[k]=%d, Typevar[k]=%d\n ",k, Tposprod[k], Typevar[k]);
+       if(FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){ /* Needs a fixed product Product of fixed dummy (<=ncovcol) covariates For a fixed product k is higher than ncovcol V3*V2 */
+       printf("MEMORY ERRORR k=%d Tvardk[k][1]=%d, Tvardk[k][2]=%d, FixedV[Tvardk[k][1]]=%d,FixedV[Tvardk[k][2]]=%d\n ",k,Tvardk[k][1],Tvardk[k][2],FixedV[Tvardk[k][1]],FixedV[Tvardk[k][2]]);
+         Fixed[k]= 0;
+         Dummy[k]= 0;
+         ncoveff++;
+         ncovf++;
+         /* ncovv++; */
+         /* TvarVV[ncovv]=Tvardk[k][1]; */
+         /* FixedV[ncovcolt+ncovv]=0; /\* or FixedV[TvarVV[ncovv]]=0 HERE *\/ */
+         /* ncovv++; */
+         /* TvarVV[ncovv]=Tvardk[k][2]; */
+         /* FixedV[ncovcolt+ncovv]=0; /\* or FixedV[TvarVV[ncovv]]=0 HERE *\/ */
+         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 */
+       }else{/* product varying Vn * Vm without age, V1+V3+age*V1+age*V3+V1*V3 looking at V1*V3, Typevar={0, 0, 1, 1, 2}, k=5, V1 is fixed, V3 is timevary, V5 is a product  */
+         /*#  ID           V1     V2          weight               birth   death   1st    s1      V3      V4      V5       2nd  s2 */
+         /* model V1+V3+age*V1+age*V3+V1*V3 + V1*V3*age*/
+         /*  Tvar={1, 3, 1, 3, 6, 6}, the 6 comes from the fact that there are already V1, V2, V3, V4, V5 native covariates */
+         k1=Tposprod[k];  /* Position in the products of product k, Tposprod={0, 0, 0, 0, 1, 1} k1=1 first product but second time varying because of V3 */
+         ncovvt++;
+         TvarVV[ncovvt]=Tvard[k1][1];  /*  TvarVV[2]=V1 (because TvarVV[1] was V3, first time varying covariates */
+         TvarVVind[ncovvt]=k;  /*  TvarVVind[2]=5 (because TvarVVind[2] was V1*V3 at position 5 */
+         ncovvt++;
+         TvarVV[ncovvt]=Tvard[k1][2];  /*  TvarVV[3]=V3 */
+         TvarVVind[ncovvt]=k;  /*  TvarVVind[2]=5 (because TvarVVind[2] was V1*V3 at position 5 */
+         /** Fixed[k] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */
+         /** Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */
+         if(Tvard[k1][1] <=ncovcol){ /* Vn is dummy fixed, (Tvard[1][1]=V1), (Tvard[1][1]=V3 time varying) */
+           if(Tvard[k1][2] <=ncovcol){ /* Vm is dummy fixed */
+             Fixed[k]= 1;
+             Dummy[k]= 0;
+             modell[k].maintype= FTYPE;
+             modell[k].subtype= FPDD;            /*      Product fixed dummy * fixed dummy */
+             ncovf++; /* Fixed variables without age */
+             TvarF[ncovf]=Tvar[k];
+             TvarFind[ncovf]=k;
+           }else if(Tvard[k1][2] <=ncovcol+nqv){ /* Vm is quanti fixed */
+             Fixed[k]= 0;  /* Fixed product */
+             Dummy[k]= 1;
+             modell[k].maintype= FTYPE;
+             modell[k].subtype= FPDQ;            /*      Product fixed dummy * fixed quantitative */
+             ncovf++; /* Varying variables without age */
+             TvarF[ncovf]=Tvar[k];
+             TvarFind[ncovf]=k;
+           }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ /* Vm is a time varying dummy covariate */
+             Fixed[k]= 1;
+             Dummy[k]= 0;
+             modell[k].maintype= VTYPE;
+             modell[k].subtype= VPDD;            /*      Product fixed dummy * varying dummy */
+             ncovv++; /* Varying variables without age */
+             TvarV[ncovv]=Tvar[k];  /* TvarV[1]=Tvar[5]=5 because there is a V4 */
+             TvarVind[ncovv]=k;/* TvarVind[1]=5 */
+           }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ /* Vm is a time varying quantitative covariate */
+             Fixed[k]= 1;
+             Dummy[k]= 1;
+             modell[k].maintype= VTYPE;
+             modell[k].subtype= VPDQ;            /*      Product fixed dummy * varying quantitative */
+             ncovv++; /* Varying variables without age */
+             TvarV[ncovv]=Tvar[k];
+             TvarVind[ncovv]=k;
+           }
+         }else if(Tvard[k1][1] <=ncovcol+nqv){ /* Vn is fixed quanti  */
+           if(Tvard[k1][2] <=ncovcol){ /* Vm is fixed dummy */
+             Fixed[k]= 0;  /*  Fixed product */
+             Dummy[k]= 1;
+             modell[k].maintype= FTYPE;
+             modell[k].subtype= FPDQ;            /*      Product fixed quantitative * fixed dummy */
+             ncovf++; /* Fixed variables without age */
+             TvarF[ncovf]=Tvar[k];
+             TvarFind[ncovf]=k;
+           }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ /* Vm is time varying */
+             Fixed[k]= 1;
+             Dummy[k]= 1;
+             modell[k].maintype= VTYPE;
+             modell[k].subtype= VPDQ;            /*      Product fixed quantitative * varying dummy */
+             ncovv++; /* Varying variables without age */
+             TvarV[ncovv]=Tvar[k];
+             TvarVind[ncovv]=k;
+           }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ /* Vm is time varying quanti */
+             Fixed[k]= 1;
+             Dummy[k]= 1;
+             modell[k].maintype= VTYPE;
+             modell[k].subtype= VPQQ;            /*      Product fixed quantitative * varying quantitative */
+             ncovv++; /* Varying variables without age */
+             TvarV[ncovv]=Tvar[k];
+             TvarVind[ncovv]=k;
+             ncovv++; /* Varying variables without age */
+             TvarV[ncovv]=Tvar[k];
+             TvarVind[ncovv]=k;
+           }
+         }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){ /* Vn is time varying dummy */
+           if(Tvard[k1][2] <=ncovcol){
+             Fixed[k]= 1;
+             Dummy[k]= 1;
+             modell[k].maintype= VTYPE;
+             modell[k].subtype= VPDD;            /*      Product time varying dummy * fixed dummy */
+             ncovv++; /* Varying variables without age */
+             TvarV[ncovv]=Tvar[k];
+             TvarVind[ncovv]=k;
+           }else if(Tvard[k1][2] <=ncovcol+nqv){
+             Fixed[k]= 1;
+             Dummy[k]= 1;
+             modell[k].maintype= VTYPE;
+             modell[k].subtype= VPDQ;            /*      Product time varying dummy * fixed quantitative */
+             ncovv++; /* Varying variables without age */
+             TvarV[ncovv]=Tvar[k];
+             TvarVind[ncovv]=k;
+           }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+             Fixed[k]= 1;
+             Dummy[k]= 0;
+             modell[k].maintype= VTYPE;
+             modell[k].subtype= VPDD;            /*      Product time varying dummy * time varying dummy */
+             ncovv++; /* Varying variables without age */
+             TvarV[ncovv]=Tvar[k];
+             TvarVind[ncovv]=k;
+           }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+             Fixed[k]= 1;
+             Dummy[k]= 1;
+             modell[k].maintype= VTYPE;
+             modell[k].subtype= VPDQ;            /*      Product time varying dummy * time varying quantitative */
+             ncovv++; /* Varying variables without age */
+             TvarV[ncovv]=Tvar[k];
+             TvarVind[ncovv]=k;
+           }
+         }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){ /* Vn is time varying quanti */
+           if(Tvard[k1][2] <=ncovcol){
+             Fixed[k]= 1;
+             Dummy[k]= 1;
+             modell[k].maintype= VTYPE;
+             modell[k].subtype= VPDQ;            /*      Product time varying quantitative * fixed dummy */
+             ncovv++; /* Varying variables without age */
+             TvarV[ncovv]=Tvar[k];
+             TvarVind[ncovv]=k;
+           }else if(Tvard[k1][2] <=ncovcol+nqv){
+             Fixed[k]= 1;
+             Dummy[k]= 1;
+             modell[k].maintype= VTYPE;
+             modell[k].subtype= VPQQ;            /*      Product time varying quantitative * fixed quantitative */
+             ncovv++; /* Varying variables without age */
+             TvarV[ncovv]=Tvar[k];
+             TvarVind[ncovv]=k;
+           }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+             Fixed[k]= 1;
+             Dummy[k]= 1;
+             modell[k].maintype= VTYPE;
+             modell[k].subtype= VPDQ;            /*      Product time varying quantitative * time varying dummy */
+             ncovv++; /* Varying variables without age */
+             TvarV[ncovv]=Tvar[k];
+             TvarVind[ncovv]=k;
+           }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+             Fixed[k]= 1;
+             Dummy[k]= 1;
+             modell[k].maintype= VTYPE;
+             modell[k].subtype= VPQQ;            /*      Product time varying quantitative * time varying quantitative */
+             ncovv++; /* Varying variables without age */
+             TvarV[ncovv]=Tvar[k];
+             TvarVind[ncovv]=k;
+           }
+         }else{
+           printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
+           fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
+         } /*end k1*/
+       }
+     }else if(Typevar[k] == 3){  /* product Vn * Vm with age, V1+V3+age*V1+age*V3+V1*V3 looking at V1*V3, Typevar={0, 0, 1, 1, 2}, k=5, V1 is fixed, V3 is timevary, V5 is a product  */
        /*#  ID           V1     V2          weight               birth   death   1st    s1      V3      V4      V5       2nd  s2 */
-       /* model V1+V3+age*V1+age*V3+V1*V3 */
+       /* model V1+V3+age*V1+age*V3+V1*V3 + V1*V3*age*/
-       /*  Tvar={1, 3, 1, 3, 6}, the 6 comes from the fact that there are already V1, V2, V3, V4, V5 native covariates */
+       /*  Tvar={1, 3, 1, 3, 6, 6}, the 6 comes from the fact that there are already V1, V2, V3, V4, V5 native covariates */
-       k1=Tposprod[k];  /* Position in the products of product k, Tposprod={0, 0, 0, 0, 1} k1=1 first product but second time varying because of V3 */
+       k1=Tposprod[k];  /* Position in the products of product k, Tposprod={0, 0, 0, 0, 1, 1} k1=1 first product but second time varying because of V3 */
-       ncovvt++;
+       ncova++;
-       TvarVV[ncovvt]=Tvard[k1][1];  /*  TvarVV[2]=V1 (because TvarVV[1] was V3, first time varying covariates */
+       TvarA[ncova]=Tvard[k1][1];  /*  TvarVV[2]=V1 (because TvarVV[1] was V3, first time varying covariates */
-       TvarVVind[ncovvt]=k;  /*  TvarVVind[2]=5 (because TvarVVind[2] was V1*V3 at position 5 */
+       TvarAind[ncova]=k;  /*  TvarVVind[2]=5 (because TvarVVind[2] was V1*V3 at position 5 */
-       ncovvt++;
+       ncova++;
-       TvarVV[ncovvt]=Tvard[k1][2];  /*  TvarVV[3]=V3 */
+       TvarA[ncova]=Tvard[k1][2];  /*  TvarVV[3]=V3 */
-       TvarVVind[ncovvt]=k;  /*  TvarVVind[2]=5 (because TvarVVind[2] was V1*V3 at position 5 */
+       TvarAind[ncova]=k;  /*  TvarVVind[2]=5 (because TvarVVind[2] was V1*V3 at position 5 */
+       /** Fixed[k] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */
+       /** Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */
+       if( FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){
+         ncovta++;
+         TvarAVVA[ncovta]=Tvard[k1][1]; /*   age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 */
+         TvarAVVAind[ncovta]=k;
+         ncovta++;
+         TvarAVVA[ncovta]=Tvard[k1][2]; /*   age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 */
+         TvarAVVAind[ncovta]=k;
+       }else{
+         ncovva++;  /* HERY  reached */
+         TvarVVA[ncovva]=Tvard[k1][1]; /*  age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4  */
+         TvarVVAind[ncovva]=k;
+         ncovva++;
+         TvarVVA[ncovva]=Tvard[k1][2]; /*   */
+         TvarVVAind[ncovva]=k;
+         ncovta++;
+         TvarAVVA[ncovta]=Tvard[k1][1]; /*   age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 */
+         TvarAVVAind[ncovta]=k;
+         ncovta++;
+         TvarAVVA[ncovta]=Tvard[k1][2]; /*   age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 */
+         TvarAVVAind[ncovta]=k;
+       }
        if(Tvard[k1][1] <=ncovcol){ /* Vn is dummy fixed, (Tvard[1][1]=V1), (Tvard[1][1]=V3 time varying) */
          if(Tvard[k1][2] <=ncovcol){ /* Vm is dummy fixed */
-           Fixed[k]= 1;
+           Fixed[k]= 2;
-           Dummy[k]= 0;
+           Dummy[k]= 2;
            modell[k].maintype= FTYPE;
            modell[k].subtype= FPDD;              /*      Product fixed dummy * fixed dummy */
-           ncovf++; /* Fixed variables without age */
+           /* TvarF[ncova]=Tvar[k];   /\* Problem to solve *\/ */
-           TvarF[ncovf]=Tvar[k];
+           /* TvarFind[ncova]=k; */
-           TvarFind[ncovf]=k;
          }else if(Tvard[k1][2] <=ncovcol+nqv){ /* Vm is quanti fixed */
-           Fixed[k]= 0;  /* Fixed product */
+           Fixed[k]= 2;  /* Fixed product */
-           Dummy[k]= 1;
+           Dummy[k]= 3;
            modell[k].maintype= FTYPE;
            modell[k].subtype= FPDQ;              /*      Product fixed dummy * fixed quantitative */
-           ncovf++; /* Varying variables without age */
+           /* TvarF[ncova]=Tvar[k]; */
-           TvarF[ncovf]=Tvar[k];
+           /* TvarFind[ncova]=k; */
-           TvarFind[ncovf]=k;
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ /* Vm is a time varying dummy covariate */
-           Fixed[k]= 1;
+           Fixed[k]= 3;
-           Dummy[k]= 0;
+           Dummy[k]= 2;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDD;              /*      Product fixed dummy * varying dummy */
-           ncovv++; /* Varying variables without age */
+           TvarV[ncova]=Tvar[k];  /* TvarV[1]=Tvar[5]=5 because there is a V4 */
-           TvarV[ncovv]=Tvar[k];  /* TvarV[1]=Tvar[5]=5 because there is a V4 */
+           TvarVind[ncova]=k;/* TvarVind[1]=5 */
-           TvarVind[ncovv]=k;/* TvarVind[1]=5 */
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ /* Vm is a time varying quantitative covariate */
-           Fixed[k]= 1;
+           Fixed[k]= 3;
-           Dummy[k]= 1;
+           Dummy[k]= 3;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDQ;              /*      Product fixed dummy * varying quantitative */
-           ncovv++; /* Varying variables without age */
+           /* ncovv++; /\* Varying variables without age *\/ */
-           TvarV[ncovv]=Tvar[k];
+           /* TvarV[ncovv]=Tvar[k]; */
-           TvarVind[ncovv]=k;
+           /* TvarVind[ncovv]=k; */
          }
        }else if(Tvard[k1][1] <=ncovcol+nqv){ /* Vn is fixed quanti  */
          if(Tvard[k1][2] <=ncovcol){ /* Vm is fixed dummy */
-           Fixed[k]= 0;  /*  Fixed product */
+           Fixed[k]= 2;  /*  Fixed product */
-           Dummy[k]= 1;
+           Dummy[k]= 2;
            modell[k].maintype= FTYPE;
            modell[k].subtype= FPDQ;              /*      Product fixed quantitative * fixed dummy */
-           ncovf++; /* Fixed variables without age */
+           /* ncova++; /\* Fixed variables with age *\/ */
-           TvarF[ncovf]=Tvar[k];
+           /* TvarF[ncovf]=Tvar[k]; */
-           TvarFind[ncovf]=k;
+           /* TvarFind[ncovf]=k; */
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ /* Vm is time varying */
-           Fixed[k]= 1;
+           Fixed[k]= 2;
-           Dummy[k]= 1;
+           Dummy[k]= 3;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDQ;              /*      Product fixed quantitative * varying dummy */
-           ncovv++; /* Varying variables without age */
+           /* ncova++; /\* Varying variables with age *\/ */
-           TvarV[ncovv]=Tvar[k];
+           /* TvarV[ncova]=Tvar[k]; */
-           TvarVind[ncovv]=k;
+           /* TvarVind[ncova]=k; */
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ /* Vm is time varying quanti */
-           Fixed[k]= 1;
+           Fixed[k]= 3;
-           Dummy[k]= 1;
+           Dummy[k]= 2;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPQQ;              /*      Product fixed quantitative * varying quantitative */
-           ncovv++; /* Varying variables without age */
+           ncova++; /* Varying variables without age */
-           TvarV[ncovv]=Tvar[k];
+           TvarV[ncova]=Tvar[k];
-           TvarVind[ncovv]=k;
+           TvarVind[ncova]=k;
-           ncovv++; /* Varying variables without age */
+           /* ncova++; /\* Varying variables without age *\/ */
-           TvarV[ncovv]=Tvar[k];
+           /* TvarV[ncova]=Tvar[k]; */
-           TvarVind[ncovv]=k;
+           /* TvarVind[ncova]=k; */
          }
        }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){ /* Vn is time varying dummy */
          if(Tvard[k1][2] <=ncovcol){
-           Fixed[k]= 1;
+           Fixed[k]= 2;
-           Dummy[k]= 1;
+           Dummy[k]= 2;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDD;              /*      Product time varying dummy * fixed dummy */
-           ncovv++; /* Varying variables without age */
+           /* ncova++; /\* Varying variables with age *\/ */
-           TvarV[ncovv]=Tvar[k];
+           /* TvarV[ncova]=Tvar[k]; */
-           TvarVind[ncovv]=k;
+           /* TvarVind[ncova]=k; */
          }else if(Tvard[k1][2] <=ncovcol+nqv){
-           Fixed[k]= 1;
+           Fixed[k]= 2;
-           Dummy[k]= 1;
+           Dummy[k]= 3;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDQ;              /*      Product time varying dummy * fixed quantitative */
-           ncovv++; /* Varying variables without age */
+           /* ncova++; /\* Varying variables with age *\/ */
-           TvarV[ncovv]=Tvar[k];
+           /* TvarV[ncova]=Tvar[k]; */
-           TvarVind[ncovv]=k;
+           /* TvarVind[ncova]=k; */
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
-           Fixed[k]= 1;
+           Fixed[k]= 3;
-           Dummy[k]= 0;
+           Dummy[k]= 2;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDD;              /*      Product time varying dummy * time varying dummy */
-           ncovv++; /* Varying variables without age */
+           /* ncova++; /\* Varying variables with age *\/ */
-           TvarV[ncovv]=Tvar[k];
+           /* TvarV[ncova]=Tvar[k]; */
-           TvarVind[ncovv]=k;
+           /* TvarVind[ncova]=k; */
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
-           Fixed[k]= 1;
+           Fixed[k]= 3;
-           Dummy[k]= 1;
+           Dummy[k]= 3;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDQ;              /*      Product time varying dummy * time varying quantitative */
-           ncovv++; /* Varying variables without age */
+           /* ncova++; /\* Varying variables with age *\/ */
-           TvarV[ncovv]=Tvar[k];
+           /* TvarV[ncova]=Tvar[k]; */
-           TvarVind[ncovv]=k;
+           /* TvarVind[ncova]=k; */
          }
        }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){ /* Vn is time varying quanti */
          if(Tvard[k1][2] <=ncovcol){
-           Fixed[k]= 1;
+           Fixed[k]= 2;
-           Dummy[k]= 1;
+           Dummy[k]= 2;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDQ;              /*      Product time varying quantitative * fixed dummy */
-           ncovv++; /* Varying variables without age */
+           /* ncova++; /\* Varying variables with age *\/ */
-           TvarV[ncovv]=Tvar[k];
+           /* TvarV[ncova]=Tvar[k]; */
-           TvarVind[ncovv]=k;
+           /* TvarVind[ncova]=k; */
          }else if(Tvard[k1][2] <=ncovcol+nqv){
-           Fixed[k]= 1;
+           Fixed[k]= 2;
-           Dummy[k]= 1;
+           Dummy[k]= 3;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPQQ;              /*      Product time varying quantitative * fixed quantitative */
-           ncovv++; /* Varying variables without age */
+           /* ncova++; /\* Varying variables with age *\/ */
-           TvarV[ncovv]=Tvar[k];
+           /* TvarV[ncova]=Tvar[k]; */
-           TvarVind[ncovv]=k;
+           /* TvarVind[ncova]=k; */
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
-           Fixed[k]= 1;
+           Fixed[k]= 3;
-           Dummy[k]= 1;
+           Dummy[k]= 2;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPDQ;              /*      Product time varying quantitative * time varying dummy */
-           ncovv++; /* Varying variables without age */
+           /* ncova++; /\* Varying variables with age *\/ */
-           TvarV[ncovv]=Tvar[k];
+           /* TvarV[ncova]=Tvar[k]; */
-           TvarVind[ncovv]=k;
+           /* TvarVind[ncova]=k; */
          }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
-           Fixed[k]= 1;
+           Fixed[k]= 3;
-           Dummy[k]= 1;
+           Dummy[k]= 3;
            modell[k].maintype= VTYPE;
            modell[k].subtype= VPQQ;              /*      Product time varying quantitative * time varying quantitative */
-           ncovv++; /* Varying variables without age */
+           /* ncova++; /\* Varying variables with age *\/ */
-           TvarV[ncovv]=Tvar[k];
+           /* TvarV[ncova]=Tvar[k]; */
-           TvarVind[ncovv]=k;
+           /* TvarVind[ncova]=k; */
          }
        }else{
          printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
          fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
        } /*end k1*/
-     }else{
+     } else{
        printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
        fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
      }
-     printf("Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[k],Dummy[k]);
+     /* printf("Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[k],Dummy[k]); */
-     printf("           modell[%d].maintype=%d, modell[%d].subtype=%d\n",k,modell[k].maintype,k,modell[k].subtype);
+     /* printf("           modell[%d].maintype=%d, modell[%d].subtype=%d\n",k,modell[k].maintype,k,modell[k].subtype); */
      fprintf(ficlog,"Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[k],Dummy[k]);
    }
+   ncovvta=ncovva;
    /* Searching for doublons in the model */
    for(k1=1; k1<= cptcovt;k1++){
      for(k2=1; k2 <k1;k2++){
- Line 11353  Dummy[k] 0=dummy (0 1), 1 quantitative (
+ Line 12254  Dummy[k] 0=dummy (0 1), 1 quantitative (
    fprintf(ficlog,"ncoveff=%d, nqfveff=%d, ntveff=%d, nqtveff=%d, cptcovn=%d\n",ncoveff,nqfveff,ntveff,nqtveff,cptcovn);
    printf("ncovf=%d, ncovv=%d, ncova=%d, nsd=%d, nsq=%d\n",ncovf,ncovv,ncova,nsd,nsq);
    fprintf(ficlog,"ncovf=%d, ncovv=%d, ncova=%d, nsd=%d, nsq=%d\n",ncovf,ncovv,ncova,nsd, nsq);
+   free_imatrix(existcomb,1,NCOVMAX,1,NCOVMAX);
    return (0); /* with covar[new additional covariate if product] and Tage if age */
    /*endread:*/
    printf("Exiting decodemodel: ");
- Line 12102  int main(int argc, char *argv[])
+ Line 13005  int main(int argc, char *argv[])
    /* double ***mobaverage; */
    double wald;
-   char line[MAXLINE];
+   char line[MAXLINE], linetmp[MAXLINE];
    char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
    char  modeltemp[MAXLINE];
- Line 12435  int main(int argc, char *argv[])
+ Line 13338  int main(int argc, char *argv[])
      }else
        break;
    }
-   if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){
+   if((num_filled=sscanf(line,"model=%[^.\n]", model)) !=EOF){ /* Every character after model but dot and  return */
+     if (num_filled != 1){
+       printf("ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);
+       fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);
+       model[0]='\0';
+       goto end;
+     }else{
+       trimbtab(linetmp,line); /* Trims multiple blanks in line */
+       strcpy(line, linetmp);
+     }
+   }
+   if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){ /* Every character after 1+age but dot and  return */
      if (num_filled != 1){
        printf("ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);
        fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);
- Line 12477  int main(int argc, char *argv[])
+ Line 13391  int main(int argc, char *argv[])
      numlinepar++;
      if(line[1]=='q'){ /* This #q will quit imach (the answer is q) */
        z[0]=line[1];
+     }else if(line[1]=='d'){ /* For debugging individual values of covariates in ficresilk */
+       debugILK=1;printf("DebugILK\n");
      }
      /* printf("****line [1] = %c \n",line[1]); */
      fputs(line, stdout);
- Line 12490  int main(int argc, char *argv[])
+ Line 13406  int main(int argc, char *argv[])
    covar=matrix(0,NCOVMAX,firstobs,lastobs);  /**< used in readdata */
    if(nqv>=1)coqvar=matrix(1,nqv,firstobs,lastobs);  /**< Fixed quantitative covariate */
    if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,firstobs,lastobs);  /**< Time varying quantitative covariate */
-   if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,firstobs,lastobs);  /**< Time varying covariate (dummy and quantitative)*/
+   /* if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,firstobs,lastobs);  /\**< Time varying covariate (dummy and quantitative)*\/ */
+   if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,ncovcol+nqv+1,ncovcol+nqv+ntv+nqtv,firstobs,lastobs);  /**< Might be better */
    cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/
    /* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5
       v1+v2*age+v2*v3 makes cptcovn = 3
- Line 12773  Please run with mle=-1 to get a correct
+ Line 13690  Please run with mle=-1 to get a correct
    TvarVQind=ivector(1,NCOVMAX); /*  */
    TvarVV=ivector(1,NCOVMAX); /*  */
    TvarVVind=ivector(1,NCOVMAX); /*  */
+   TvarVVA=ivector(1,NCOVMAX); /*  */
+   TvarVVAind=ivector(1,NCOVMAX); /*  */
+   TvarAVVA=ivector(1,NCOVMAX); /*  */
+   TvarAVVAind=ivector(1,NCOVMAX); /*  */
    Tvalsel=vector(1,NCOVMAX); /*  */
    Tvarsel=ivector(1,NCOVMAX); /*  */
    Typevar=ivector(-1,NCOVMAX); /* -1 to 2 */
    Fixed=ivector(-1,NCOVMAX); /* -1 to 3 */
    Dummy=ivector(-1,NCOVMAX); /* -1 to 3 */
+   DummyV=ivector(-1,NCOVMAX); /* 1 to 3 */
+   FixedV=ivector(-1,NCOVMAX); /* 1 to 3 */
    /*  V2+V1+V4+age*V3 is a model with 4 covariates (3 plus signs).
        For each model-covariate stores the data-covariate id. Tvar[1]=2, Tvar[2]=1, Tvar[3]=4,
        Tvar[4=age*V3] is 3 and 'age' is recorded in Tage.
- Line 12798  Please run with mle=-1 to get a correct
+ Line 13722  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);
+   Tvardk=imatrix(0,NCOVMAX,1,2);
    Tage=ivector(1,NCOVMAX); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age
 covariates (3 plus signs)
                           Tage[1=V3*age]= 4; Tage[2=age*V4] = 3
- Line 12818  Please run with mle=-1 to get a correct
+ Line 13742  Please run with mle=-1 to get a correct
                                  * Tmodelqind[1]=1,Tvaraff[1]@9={4,
                                  * 3, 1, 0, 0, 0, 0, 0, 0},
                                  * model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/
+ /* Probably useless zeroes */
+   for(i=1;i<NCOVMAX;i++){
+     DummyV[i]=0;
+     FixedV[i]=0;
+   }
+   for(i=1; i <=ncovcol;i++){
+     DummyV[i]=0;
+     FixedV[i]=0;
+   }
+   for(i=ncovcol+1; i <=ncovcol+nqv;i++){
+     DummyV[i]=1;
+     FixedV[i]=0;
+   }
+   for(i=ncovcol+nqv+1; i <=ncovcol+nqv+ntv;i++){
+     DummyV[i]=0;
+     FixedV[i]=1;
+   }
+   for(i=ncovcol+nqv+ntv+1; i <=ncovcol+nqv+ntv+nqtv;i++){
+     DummyV[i]=1;
+     FixedV[i]=1;
+   }
+   for(i=1; i <=ncovcol+nqv+ntv+nqtv;i++){
+     printf("Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",i,i,DummyV[i],i,FixedV[i]);
+     fprintf(ficlog,"Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",i,i,DummyV[i],i,FixedV[i]);
+   }
  /* Main decodemodel */
- Line 13394  Please run with mle=-1 to get a correct
+ Line 14348  Please run with mle=-1 to get a correct
          fprintf(ficres,"  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
          fprintf(ficlog,"  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
          fprintf(fichtm, "<th>+  V%d*V%d</th>",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+       }else if(Typevar[j]==3) { /* TO VERIFY */
+         printf("  +    V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+         fprintf(ficres,"  +    V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+         fprintf(ficlog,"  +    V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+         fprintf(fichtm, "<th>+  V%d*V%d*age</th>",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
        }
      }
      printf("\n");
- Line 13453  Please run with mle=-1 to get a correct
+ Line 14412  Please run with mle=-1 to get a correct
            fprintf(fichtm, "<th>+  V%d*age</th>",Tvar[j]);
          }else if(Typevar[j]==2) {
            fprintf(fichtm, "<th>+  V%d*V%d</th>",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+         }else if(Typevar[j]==3) { /* TO VERIFY */
+           fprintf(fichtm, "<th>+  V%d*V%d*age</th>",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
          }
        }
        fprintf(fichtm, "</tr>\n");
- Line 13510  Please run with mle=-1 to get a correct
+ Line 14471  Please run with mle=-1 to get a correct
      }
      fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
-     if(mle >= 1) /* To big for the screen */
+     if(mle >= 1) /* Too big for the screen */
        printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
      fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n#   ...\n# 232 Cov(b23,a12)  Cov(b23,b12) ... Var (b23)\n");
      /* # 121 Var(a12)\n\ */
- Line 13774  Please run with mle=-1 to get a correct
+ Line 14735  Please run with mle=-1 to get a correct
        case 13:
          num_filled=sscanf(line,"result:%[^\n]\n",resultlineori);
          nresult++; /* Sum of resultlines */
-         printf("Result %d: result:%s\n",nresult, resultlineori);
+         /* printf("Result %d: result:%s\n",nresult, resultlineori); */
          /* removefirstspace(&resultlineori); */
          if(strstr(resultlineori,"v") !=0){
- Line 13783  Please run with mle=-1 to get a correct
+ Line 14744  Please run with mle=-1 to get a correct
            return 1;
          }
          trimbb(resultline, resultlineori); /* Suppressing double blank in the resultline */
-         printf("Decoderesult resultline=\"%s\" resultlineori=\"%s\"\n", resultline, resultlineori);
+         /* printf("Decoderesult resultline=\"%s\" resultlineori=\"%s\"\n", resultline, resultlineori); */
          if(nresult > MAXRESULTLINESPONE-1){
            printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINESPONE-1,nresult,rfileres);
            fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINESPONE-1,nresult,rfileres);
- Line 13866  Please run with mle=-1 to get a correct
+ Line 14827  Please run with mle=-1 to get a correct
          date2dmy(datebackf,&jbackf, &mbackf, &anbackf);
        }
-       printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,bage, fage, prevfcast, prevbcast, pathc,p, (int)anprojd-bage, (int)anbackd-fage);
+       printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,bage, fage, prevfcast, prevbcast, pathc,p, (int)anprojd-bage, (int)anbackd-fage);/* HERE valgrind Tvard*/
      }
      printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \
                   model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,prevbcast, estepm, \
- Line 14016  Please run with mle=-1 to get a correct
+ Line 14977  Please run with mle=-1 to get a correct
      pstamp(ficreseij);
-     i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
+     /* i1=pow(2,cptcoveff); /\* Number of combination of dummy covariates *\/ */
-     if (cptcovn < 1){i1=1;}
+     /* if (cptcovn < 1){i1=1;} */
-     for(nres=1; nres <= nresult; nres++) /* For each resultline */
+     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-     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; */
        fprintf(ficreseij,"\n#****** ");
        printf("\n#****** ");
-       for(j=1;j<=cptcoveff;j++) {
+       for(j=1;j<=cptcovs;j++){
-         fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+       /* for(j=1;j<=cptcoveff;j++) { */
-         printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+         /* fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
+         fprintf(ficreseij," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
+         printf(" V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
+         /* printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
        }
        for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
          printf(" V%d=%lg ",TvarsQ[j], TinvDoQresult[nres][TvarsQ[j]]); /* TvarsQ[j] gives the name of the jth quantitative (fixed or time v) */
- Line 14096  Please run with mle=-1 to get a correct
+ Line 15060  Please run with mle=-1 to get a correct
        /* */
        if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */
          continue;
-       printf("\n# model %s \n#****** Result for:", model);
+       printf("\n# model %s \n#****** Result for:", model);  /* HERE model is empty */
        fprintf(ficrest,"\n# model %s \n#****** Result for:", model);
        fprintf(ficlog,"\n# model %s \n#****** Result for:", model);
        /* It might not be a good idea to mix dummies and quantitative */
- Line 14120  Please run with mle=-1 to get a correct
+ Line 15084  Please run with mle=-1 to get a correct
          /* Tvresult[nres][j] Name of the variable at position j in this resultline */
          /* Tresult[nres][j] Value of this variable at position j could be a float if quantitative  */
  /* We give up with the combinations!! */
-         printf("\n j=%d In computing T_ Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=%d cptcovs=%d, cptcoveff=%d Fixed[modelresult[nres][j]]=%d\n", j, nres, j, Dummy[modelresult[nres][j]],nres,j,modelresult[nres][j],cptcovs, cptcoveff,Fixed[modelresult[nres][j]]);  /* end if dummy  or quanti */
+         /* if(debugILK) */
+         /*   printf("\n j=%d In computing T_ Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=%d cptcovs=%d, cptcoveff=%d Fixed[modelresult[nres][j]]=%d\n", j, nres, j, Dummy[modelresult[nres][j]],nres,j,modelresult[nres][j],cptcovs, cptcoveff,Fixed[modelresult[nres][j]]);  /\* end if dummy  or quanti *\/ */
          if(Dummy[modelresult[nres][j]]==0){/* Dummy variable of the variable in position modelresult in the model corresponding to j in resultline  */
-           printf("V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
+           /* printf("V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][j]); /\* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  *\/ */ /* TinvDoQresult[nres][Name of the variable] */
-           fprintf(ficlog,"V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
+           printf("V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); /* Output of each value for the combination TKresult[nres], ordered by the covariate values in the resultline  */
-           fprintf(ficrest,"V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
+           fprintf(ficlog,"V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
+           fprintf(ficrest,"V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline  */
            if(Fixed[modelresult[nres][j]]==0){ /* Fixed */
              printf("fixed ");fprintf(ficlog,"fixed ");fprintf(ficrest,"fixed ");
            }else{
- Line 14269  Please run with mle=-1 to get a correct
+ Line 15235  Please run with mle=-1 to get a correct
      free_vector(weight,firstobs,lastobs);
-     free_imatrix(Tvardk,1,NCOVMAX,1,2);
+     free_imatrix(Tvardk,0,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);
- Line 14294  Please run with mle=-1 to get a correct
+ Line 15260  Please run with mle=-1 to get a correct
    free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
    free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
    free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
-   if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,firstobs,lastobs);
+   /* if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,firstobs,lastobs); */
+   if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,ncovcol+nqv+1,ncovcol+nqv+ntv+nqtv,firstobs,lastobs);
    if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,firstobs,lastobs);
    if(nqv>=1)free_matrix(coqvar,1,nqv,firstobs,lastobs);
    free_matrix(covar,0,NCOVMAX,firstobs,lastobs);
- Line 14310  Please run with mle=-1 to get a correct
+ Line 15277  Please run with mle=-1 to get a correct
    free_ivector(ncodemaxwundef,1,NCOVMAX);
    free_ivector(Dummy,-1,NCOVMAX);
    free_ivector(Fixed,-1,NCOVMAX);
-   free_ivector(DummyV,1,NCOVMAX);
+   free_ivector(DummyV,-1,NCOVMAX);
-   free_ivector(FixedV,1,NCOVMAX);
+   free_ivector(FixedV,-1,NCOVMAX);
    free_ivector(Typevar,-1,NCOVMAX);
    free_ivector(Tvar,1,NCOVMAX);
    free_ivector(TvarsQ,1,NCOVMAX);
- Line 14333  Please run with mle=-1 to get a correct
+ Line 15300  Please run with mle=-1 to get a correct
    free_ivector(TvarVDind,1,NCOVMAX);
    free_ivector(TvarVQ,1,NCOVMAX);
    free_ivector(TvarVQind,1,NCOVMAX);
+   free_ivector(TvarAVVA,1,NCOVMAX);
+   free_ivector(TvarAVVAind,1,NCOVMAX);
+   free_ivector(TvarVVA,1,NCOVMAX);
+   free_ivector(TvarVVAind,1,NCOVMAX);
    free_ivector(TvarVV,1,NCOVMAX);
    free_ivector(TvarVVind,1,NCOVMAX);

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

Removed from v.1.339
changed lines
	Added in v.1.352