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.346 and 1.358

-version 1.346, 2022/09/16 13:52:36
+version 1.358, 2023/06/14 14:57:02
  Line 1
  /* $Id$
    $State$
    $Log$
+   Revision 1.358  2023/06/14 14:57:02  brouard
+   * imach.c (Module): Testing if conjugate gradient could be quicker when lot of variables POWELLORIGINCONJUGATE
+   Revision 1.357  2023/06/14 14:55:52  brouard
+   * imach.c (Module): Testing if conjugate gradient could be quicker when lot of variables POWELLORIGINCONJUGATE
+   Revision 1.356  2023/05/23 12:08:43  brouard
+   Summary: 0.99r46
+   * imach.c (Module): Fixed PROB_r
+   Revision 1.355  2023/05/22 17:03:18  brouard
+   Summary: 0.99r46
+   * imach.c (Module): In the ILK....txt file, the number of columns
+   before the covariates values is dependent of the number of states (16+nlstate): 0.99r46
+   Revision 1.354  2023/05/21 05:05:17  brouard
+   Summary: Temporary change for imachprax
+   Revision 1.353  2023/05/08 18:48:22  brouard
+   *** empty log message ***
+   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
- Line 1260  Important routines
+ Line 1298  Important routines
  /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
  /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
  /* #define FLATSUP  *//* Suppresses directions where likelihood is flat */
+ #define POWELLORIGINCONJUGATE  /* Don't use conjugate but biggest decrease if valuable */
  #include <math.h>
  #include <stdio.h>
- Line 1311  typedef struct {
+ Line 1350  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
- Line 1322  double gnuplotversion=GNUPLOTVERSION;
+ Line 1361  double gnuplotversion=GNUPLOTVERSION;
  #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 1356  double gnuplotversion=GNUPLOTVERSION;
+ Line 1395  double gnuplotversion=GNUPLOTVERSION;
  /* $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[]="Testing conjugate 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];
- Line 1370  int cptcovt=0; /**< cptcovt: total numbe
+ Line 1409  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 1462  extern time_t time();
+ Line 1507  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 1565  int **nbcode, *Tvar; /**< model=V2 => Tv
+ Line 1611  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 */
- /*    kmodel  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 1635  int *TvarVQ; /* TvarVQ[1]=V5 in V5+V4+V3
+ Line 1686  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 1779  char *trimbb(char *out, char *in)
+ Line 1834  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 1805  char *trimbb(char *out, char *in)
+ Line 1874  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 1825  char *substrchaine(char *out, char *in,
+ Line 1894  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 1891  int nbocc(char *s, char occ)
+ Line 1960  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 2570  void powell(double p[], double **xi, int
+ Line 2661  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 2583  void powell(double p[], double **xi, int
+ Line 2675  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 2612  void powell(double p[], double **xi, int
+ Line 2706  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 2642  void powell(double p[], double **xi, int
+ Line 2739  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 2663  void powell(double p[], double **xi, int
+ Line 2762  void powell(double p[], double **xi, int
        printf("%d",i);fflush(stdout); /* print direction (parameter) i */
        fprintf(ficlog,"%d",i);fflush(ficlog);
  #ifdef LINMINORIGINAL
-       linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input.*/
+       linmin(p,xit,n,fret,func); /* New point i minimizing in direction i has coordinates p[j].*/
+       /* xit[j] gives the n coordinates of direction i as input.*/
+       /* *fret gives the maximum value on direction xit */
  #else
        linmin(p,xit,n,fret,func,&flat); /* Point p[n]. xit[n] has been loaded for direction i as input.*/
                          flatdir[i]=flat; /* Function is vanishing in that direction i */
- Line 2693  void powell(double p[], double **xi, int
+ Line 2794  void powell(double p[], double **xi, int
        fprintf(ficlog,"\n");
  #endif
      } /* end loop on each direction i */
-     /* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */
+     /* Convergence test will use last linmin estimation (fret) and compare to former iteration (fp) */
      /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit  */
-     /* New value of last point Pn is not computed, P(n-1) */
      for(j=1;j<=n;j++) {
        if(flatdir[j] >0){
          printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
- Line 2833  void powell(double p[], double **xi, int
+ Line 2933  void powell(double p[], double **xi, int
            }
          }
  #endif
+ #ifdef POWELLORIGINCONJUGATE
          for (j=1;j<=n;j++) {
            xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
            xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
          }
+ #else
+         for (j=1;j<=n-1;j++) {
+           xi[j][1]=xi[j][j+1]; /* Standard method of conjugate directions */
+           xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
+         }
+ #endif
  #ifdef LINMINORIGINAL
  #else
          for (j=1, flatd=0;j<=n;j++) {
- Line 2865  void powell(double p[], double **xi, int
+ Line 2972  void powell(double p[], double **xi, int
  #endif
          printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
          fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
+   /* The minimization in direction $\xi_1$ gives $P_1$. From $P_1$ minimization in direction $\xi_2$ gives */
+   /* $P_2$. Minimization of line $P_2-P_1$ gives new starting point $P^{(1)}_0$ and direction $\xi_1$ is dropped and replaced by second */
+   /* direction $\xi_1^{(1)}=\xi_2$. Also second direction is replaced by new direction $\xi^{(1)}_2=P_2-P_0$. */
+   /* At the second iteration, starting from $P_0^{(1)}$, minimization amongst $\xi^{(1)}_1$ gives point $P^{(1)}_1$. */
+   /* Minimization amongst $\xi^{(1)}_2=(P_2-P_0)$ gives point $P^{(1)}_2$.  As $P^{(2)}_1$ and */
+   /* $P^{(1)}_0$ are minimizing in the same direction $P^{(1)}_2 - P^{(1)}_1= P_2-P_0$, directions $P^{(1)}_2-P^{(1)}_0$ */
+   /* and $P_2-P_0$ (parallel to $\xi$ and $\xi^c$) are conjugate.  } */
  #ifdef DEBUG
          printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
- Line 2951  void powell(double p[], double **xi, int
+ Line 3067  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 3161  void powell(double p[], double **xi, int
+ Line 3277  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 3628  double ***hpxij(double ***po, int nhstep
+ Line 3744  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 3814  double ***hbxij(double ***po, int nhstep
+ Line 3930  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 3999  double func( double *x)
+ Line 4115  double func( double *x)
            iposold=ipos;
            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 4026  double func( double *x)
+ Line 4130  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; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
+           /*     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 4300  double func( double *x)
+ Line 4422  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;
- Line 4336  double funcone( double *x)
+ Line 4458  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 4393  double funcone( double *x)
+ Line 4515  double funcone( double *x)
        * 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[4]=6
+        * 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
-         *             V1  V2     V3    V4   V5 V6     V7  V8
+         * 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      0      0    0  1      1   1  0, 0, 1,1,   1, 0, 1, 0, 1, 0, 1, 0}
-         * model=          V2  +  V3  +  V4  +  V6  +  V7  +  V6*V2  +  V7*V2  +  V6*V3  +  V7*V3  +  V6*V4  +  V7*V4
+         * 3           0   0      0      0    0  1      1   1  0,     1      1    1      1    1}
-         * kmodel           1     2      3      4      5        6         7         8         9        10        11
+         * 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
-         * ncovvt=14                            1      2       3 4       5 6       7 8       9 10     11 12     13 14
+         *                                    V6       V7      V6*V2     V7*V2     V6*V3     V7*V3     V6*V4     V7*V4
-         * TvarVV[1]@14 = itv                   {6,     7,     6, 2,     7, 2,     6, 3,     7, 3,     6, 4,     7, 4}
+         * ncovvt=14                            1      2        3 4       5 6       7 8       9 10     11 12     13 14
-         * TvarVVind[1]@14=                    {4,     5,      6, 6,     7, 7,     8, 8,      9, 9,   10, 10,   11, 11}
+         * 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]@20=     {2,    3,     4,    6,      7,      9,      10,        11,       12,      13,       14}
+         * 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
+         *? FixedV[itv]                          1      1       1  0      1 0       1 0       1 0      1 0     1 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++){ /* Varying  covariates (single and product but no age) including individual from products */
+       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, exploding product  */
+         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 */
+           /* printf("DEBUG ncovv=%d, Varying TvarVV[ncovv]=%d\n",ncovv, TvarVV[ncovv]); */
            cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i];  /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */
+           /* printf("DEBUG Varying cov[ioffset+ipos=%d]=%g \n",ioffset+ipos,cotvarv); */
          }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];  /* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model */
+           /* printf("DEBUG ncovv=%d, Fixed TvarVV[ncovv]=%d\n",ncovv, TvarVV[ncovv]); */
+           cotvarv=covar[itv][i];  /* Good: In V6*V3, 3 is fixed at position of the data */
+           /* printf("DEBUG Fixed cov[ioffset+ipos=%d]=%g \n",ioffset+ipos,cotvarv); */
          }
          if(ipos!=iposold){ /* Not a product or first of a product */
            cotvarvold=cotvarv;
- Line 4430  double funcone( double *x)
+ Line 4634  double funcone( double *x)
          }
          iposold=ipos;
          cov[ioffset+ipos]=cotvarv;
+         /* printf("DEBUG Product cov[ioffset+ipos=%d] \n",ioffset+ipos); */
          /* For products */
        }
        /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates single *\/ */
- Line 4468  double funcone( double *x)
+ Line 4673  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; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
+           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 4561  double funcone( double *x)
+ Line 4784  double funcone( double *x)
            }
            iposold=ipos;
          }
-         for (kk=1; kk<=cptcovage;kk++) {
+         /* for (kk=1; kk<=cptcovage;kk++) { */
-           if(!FixedV[Tvar[Tage[kk]]]){
+         /*   if(!FixedV[Tvar[Tage[kk]]]){ */
-             fprintf(ficresilk," %g*age",covar[Tvar[Tage[kk]]][i]);
+         /*     fprintf(ficresilk," %g*age",covar[Tvar[Tage[kk]]][i]); */
-             /* printf(" %g*age",covar[Tvar[Tage[kk]]][i]); */
+         /*     /\* printf(" %g*age",covar[Tvar[Tage[kk]]][i]); *\/ */
-           }else{
+         /*   }else{ */
-             fprintf(ficresilk," %g*age",cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]);/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
+         /*     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) *\/  */
+         /*     /\* 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 *\/ */
- Line 4662  void likelione(FILE *ficres,double p[],
+ Line 4910  void likelione(FILE *ficres,double p[],
        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 */
+          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): <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
+          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]]);
- <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]]);
+          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 */
- Line 4728  void mlikeli(FILE *ficres,double p[], in
+ Line 4977  void mlikeli(FILE *ficres,double p[], in
    double fret;
    double fretone; /* Only one call to likelihood */
    /*  char filerespow[FILENAMELENGTH];*/
+   double * p1; /* Shifted parameters from 0 instead of 1 */
  #ifdef NLOPT
    int creturn;
    nlopt_opt opt;
    /* double lb[9] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL }; /\* lower bounds *\/ */
    double *lb;
    double minf; /* the minimum objective value, upon return */
-   double * p1; /* Shifted parameters from 0 instead of 1 */
    myfunc_data dinst, *d = &dinst;
  #endif
    xi=matrix(1,npar,1,npar);
-   for (i=1;i<=npar;i++)
+   for (i=1;i<=npar;i++)  /* Starting with canonical directions j=1,n xi[i=1,n][j] */
      for (j=1;j<=npar;j++)
        xi[i][j]=(i==j ? 1.0 : 0.0);
    printf("Powell\n");  fprintf(ficlog,"Powell\n");
- Line 5070  double hessij( double x[], double **hess
+ Line 5320  double hessij( double x[], double **hess
        kmax=kmax+10;
      }
      if(kmax >=10 || firstime ==1){
+       /* What are the thetai and thetaj? thetai/ncovmodel thetai=(thetai-thetai%ncovmodel)/ncovmodel +thetai%ncovmodel=(line,pos)  */
        printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
        fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
        printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
- Line 6284  void  concatwav(int wav[], int **dh, int
+ Line 6535  void  concatwav(int wav[], int **dh, int
     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); */
-      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 6381  void  concatwav(int wav[], int **dh, int
+ Line 6632  void  concatwav(int wav[], int **dh, int
           break;
         } /* end switch */
       } /* end dummy test */
-      if(Dummy[k]==1 && Typevar[k] !=1 && Fixed ==0){ /* Fixed 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 7407  void varprob(char optionfilefiname[], do
+ Line 7658  void varprob(char optionfilefiname[], do
     double ***varpij;
     strcpy(fileresprob,"PROB_");
-    strcat(fileresprob,fileres);
+    strcat(fileresprob,fileresu);
     if((ficresprob=fopen(fileresprob,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", fileresprob);
       fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
- Line 7564  To be simple, these graphs help to under
+ Line 7815  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 8155  true period expectancies (those weighted
+ Line 8406  true period expectancies (those weighted
  /******************* Gnuplot file **************/
  void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double bage, double fage , int prevfcast, int prevbcast, char pathc[], double p[], int offyear, int offbyear){
-   char dirfileres[132],optfileres[132];
+   char dirfileres[256],optfileres[256];
-   char gplotcondition[132], gplotlabel[132];
+   char gplotcondition[256], gplotlabel[256];
    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;
- Line 8224  void printinggnuplot(char fileresu[], ch
+ Line 8475  void printinggnuplot(char fileresu[], ch
    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 */
+     kvar=Tvar[TvarFind[kf]]; /* variable name */
-     k=18+Tvar[TvarFind[kf]];/*offset because there are 18 columns in the ILK_ file */
+     /* 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 *\/ */
+     k=16+nlstate+kf;/*offset because there are 19 columns in the ILK_ file, first cov Vn on col 21 with 4 living states */
      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));
-       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);
+       if(gnuplotversion >=5.2){ /* Former gnuplot versions do not have variable pointsize!! */
-       for (j=2; j<= nlstate+ndeath ; j ++) {
+         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);
-         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);
+         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");
      }
- Line 8261  void printinggnuplot(char fileresu[], ch
+ Line 8522  void printinggnuplot(char fileresu[], ch
            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);
- Line 8903  set ter svg size 640, 480\nunset log y\n
+ Line 9165  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 8998  set ter svg size 640, 480\nunset log y\n
+ Line 9264  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 9110  set ter svg size 640, 480\nunset log y\n
+ Line 9377  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 9235  set ter svg size 640, 480\nunset log y\n
+ Line 9503  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 9321  set ter svg size 640, 480\nunset log y\n
+ Line 9614  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 9674  void prevforecast(char fileres[], double
+ Line 9996  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 9722  void prevforecast(char fileres[], double
+ Line 10050  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 9816  void prevforecast(char fileres[], double
+ Line 10146  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 9869  void prevforecast(char fileres[], double
+ Line 10205  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 10471  int readdata(char datafile[], int firsto
+ Line 10809  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 10650  int readdata(char datafile[], int firsto
+ Line 10963  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 10855  int decoderesult( char resultline[], int
+ Line 11172  int decoderesult( char resultline[], int
    if (strlen(resultsav) >1){
      j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' in this resultline */
    }
-   if(j == 0){ /* Resultline but no = */
+   if(j == 0 && cptcovs== 0){ /* Resultline but no =  and no covariate in the model */
      TKresult[nres]=0; /* Combination for the nresult and the model */
      return (0);
    }
    if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */
-     printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model);
+     fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, 1+age+%s.\n",j, cptcovs, model);fflush(ficlog);
-     fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model);
+     printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, 1+age+%s.\n",j, cptcovs, model);fflush(stdout);
-     /* return 1;*/
+     if(j==0)
+       return 1;
    }
    for(k=1; k<=j;k++){ /* Loop on any covariate of the RESULT LINE */
      if(nbocc(resultsav,'=') >1){
- Line 10921  int decoderesult( char resultline[], int
+ Line 11239  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]); */
- Line 10933  int decoderesult( char resultline[], int
+ Line 11251  int decoderesult( char resultline[], int
          }
        }
        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;
- Line 10947  int decoderesult( char resultline[], int
+ Line 11265  int decoderesult( char resultline[], int
          }
        }
        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 10960  int decoderesult( char resultline[], int
+ Line 11278  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 11045  int decoderesult( char resultline[], int
+ Line 11363  int decoderesult( char resultline[], int
        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]); */
        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*/
+         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{
+         k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 25.1=k3q */
+         k2q=(int)Tvarsel[k3q]; /*  Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */
+         TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* 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 ){ /* For product quant or dummy (not with age) */
+     }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]]); */
      }else{
- Line 11092  int decodemodel( char model[], int lasto
+ Line 11419  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 11131  int decodemodel( char model[], int lasto
+ Line 11465  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 11197  int decodemodel( char model[], int lasto
+ Line 11533  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 11204  int decodemodel( char model[], int lasto
+ Line 11556  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;
+               }
+             }else{  /* k1 Vn*Vm already exists */
+               k11=existcomb[n][m];
+               Tposprod[k]=k11; /* OK */
+               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 */
                }
-             } /*End of FixedV */
+             }else{ /*  for product Vn*Vm */
-           } /* End age is not in the model */
+               Typevar[k]=2;  /* 2 for product Vn*Vm */
-         } /* End if model includes a product */
+               cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
-         else { /* not a product */
+               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 11276  int decodemodel( char model[], int lasto
+ Line 11758  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 11304  int decodemodel( char model[], int lasto
+ Line 11789  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<=NCOVMAX; k++){ Fixed[k]=0; Dummy[k]=0;}
    for(k=1;k<=NCOVMAX; k++){TvarFind[k]=0; TvarVind[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 loop on k from model */
+   /* 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 11329  Dummy[k] 0=dummy (0 1), 1 quantitative (
+ Line 11819  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 11401  Dummy[k] 0=dummy (0 1), 1 quantitative (
+ Line 11880  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); */
      }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]);
      }
- Line 11588  Dummy[k] 0=dummy (0 1), 1 quantitative (
+ Line 12277  Dummy[k] 0=dummy (0 1), 1 quantitative (
      /* 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 11615  Dummy[k] 0=dummy (0 1), 1 quantitative (
+ Line 12305  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 12364  int main(int argc, char *argv[])
+ Line 13056  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 12697  int main(int argc, char *argv[])
+ Line 13389  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 13038  Please run with mle=-1 to get a correct
+ Line 13741  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 13063  Please run with mle=-1 to get a correct
+ Line 13773  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 13083  Please run with mle=-1 to get a correct
+ Line 13793  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 13659  Please run with mle=-1 to get a correct
+ Line 14399  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 13718  Please run with mle=-1 to get a correct
+ Line 14463  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 13775  Please run with mle=-1 to get a correct
+ Line 14522  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 14131  Please run with mle=-1 to get a correct
+ Line 14878  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 14281  Please run with mle=-1 to get a correct
+ Line 15028  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 14361  Please run with mle=-1 to get a correct
+ Line 15111  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 14536  Please run with mle=-1 to get a correct
+ Line 15286  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 14578  Please run with mle=-1 to get a correct
+ Line 15328  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 14601  Please run with mle=-1 to get a correct
+ Line 15351  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.346
changed lines
	Added in v.1.358