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.209 and 1.218

-version 1.209, 2015/11/17 22:12:03
+version 1.218, 2016/02/12 11:29:23
  Line 1
  /* $Id$
    $State$
    $Log$
+   Revision 1.218  2016/02/12 11:29:23  brouard
+   Summary: 0.99 Back projections
+   Revision 1.217  2015/12/23 17:18:31  brouard
+   Summary: Experimental backcast
+   Revision 1.216  2015/12/18 17:32:11  brouard
+   Summary: 0.98r4 Warning and status=-2
+   Version 0.98r4 is now:
+    - displaying an error when status is -1, date of interview unknown and date of death known;
+    - permitting a status -2 when the vital status is unknown at a known date of right truncation.
+   Older changes concerning s=-2, dating from 2005 have been supersed.
+   Revision 1.215  2015/12/16 08:52:24  brouard
+   Summary: 0.98r4 working
+   Revision 1.214  2015/12/16 06:57:54  brouard
+   Summary: temporary not working
+   Revision 1.213  2015/12/11 18:22:17  brouard
+   Summary: 0.98r4
+   Revision 1.212  2015/11/21 12:47:24  brouard
+   Summary: minor typo
+   Revision 1.211  2015/11/21 12:41:11  brouard
+   Summary: 0.98r3 with some graph of projected cross-sectional
+   Author: Nicolas Brouard
+   Revision 1.210  2015/11/18 17:41:20  brouard
+   Summary: Start working on projected prevalences
    Revision 1.209  2015/11/17 22:12:03  brouard
    Summary: Adding ftolpl parameter
    Author: N Brouard
- Line 620
+ Line 654
    hPijx.
    Also this programme outputs the covariance matrix of the parameters but also
    of the life expectancies. It also computes the period (stable) prevalence.
+ Back prevalence and projections:
+  - back_prevalence_limit(double *p, double **bprlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp, double dateprev1,double dateprev2, int firstpass, int lastpass, int mobilavproj)
+     Computes the back prevalence limit  for any combination     of covariate values k
+     at any age between ageminpar and agemaxpar and returns it in **bprlim. In the loops,
+    - **bprevalim(**bprlim, ***mobaverage, nlstate, *p, age, **oldm, **savm, **dnewm, **doldm, **dsavm, ftolpl, ncvyearp, k);
+  - hBijx Back Probability to be in state i at age x-h being in j at x
+    Computes for any combination of covariates k and any age between bage and fage
+    p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+                         oldm=oldms;savm=savms;
+          - hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
+      Computes the transition matrix starting at age 'age' over
+      'nhstepm*hstepm*stepm' months (i.e. until
+      age (in years)  age+nhstepm*hstepm*stepm/12) by multiplying
+      nhstepm*hstepm matrices. Returns p3mat[i][j][h] after calling
+      p3mat[i][j][h]=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\
+,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
    Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).
             Institut national d'études démographiques, Paris.
    This software have been partly granted by Euro-REVES, a concerted action
- Line 754  typedef struct {
+ Line 805  typedef struct {
  #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
  #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */
  #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
+ /*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/
+ #define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1
  #define MAXN 20000
  #define YEARM 12. /**< Number of months per year */
- #define AGESUP 130
+ /* #define AGESUP 130 */
+ #define AGESUP 150
+ #define AGEMARGE 25 /* Marge for agemin and agemax for(iage=agemin-AGEMARGE; iage <= agemax+3+AGEMARGE; iage++) */
  #define AGEBASE 40
  #define AGEOVERFLOW 1.e20
  #define AGEGOMP 10 /**< Minimal age for Gompertz adjustment */
- Line 788  int cptcovage=0; /**< Number of covariat
+ Line 843  int cptcovage=0; /**< Number of covariat
  int cptcovprodnoage=0; /**< Number of covariate products without age */
  int cptcoveff=0; /* Total number of covariates to vary for printing results */
  int cptcov=0; /* Working variable */
+ int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */
  int npar=NPARMAX;
  int nlstate=2; /* Number of live states */
  int ndeath=1; /* Number of dead states */
- Line 810  double jmean=1; /* Mean space between 2
+ Line 866  double jmean=1; /* Mean space between 2
  double **matprod2(); /* test */
  double **oldm, **newm, **savm; /* Working pointers to matrices */
  double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
+ double   **ddnewms, **ddoldms, **ddsavms; /* for freeing later */
  /*FILE *fic ; */ /* Used in readdata only */
- FILE *ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
+ FILE *ficpar, *ficparo,*ficres, *ficresp, *ficresphtm, *ficresphtmfr, *ficrespl, *ficresplb,*ficrespij, *ficrespijb, *ficrest,*ficresf, *ficresfb,*ficrespop;
  FILE *ficlog, *ficrespow;
  int globpr=0; /* Global variable for printing or not */
  double fretone; /* Only one call to likelihood */
- Line 834  char fileresv[FILENAMELENGTH];
+ Line 892  char fileresv[FILENAMELENGTH];
  FILE  *ficresvpl;
  char fileresvpl[FILENAMELENGTH];
  char title[MAXLINE];
- char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];
+ char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH],  fileresplb[FILENAMELENGTH];
  char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH];
  char tmpout[FILENAMELENGTH],  tmpout2[FILENAMELENGTH];
  char command[FILENAMELENGTH];
  int  outcmd=0;
- char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
+ char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filerespijb[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
  char fileresu[FILENAMELENGTH]; /* fileres without r in front */
  char filelog[FILENAMELENGTH]; /* Log file */
  char filerest[FILENAMELENGTH];
- Line 917  int *ncodemaxwundef;  /* ncodemax[j]= Nu
+ Line 975  int *ncodemaxwundef;  /* ncodemax[j]= Nu
                               covariate for which somebody answered including
                               undefined. Usually 3: -1, 0 and 1. */
  double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
- double **pmmij, ***probs;
+ double **pmmij, ***probs; /* Global pointer */
+ double ***mobaverage; /* New global variable */
  double *ageexmed,*agecens;
  double dateintmean=0;
- Line 1373  char *subdirf3(char fileres[], char *pre
+ Line 1432  char *subdirf3(char fileres[], char *pre
    strcat(tmpout,fileres);
    return tmpout;
  }
+ /*************** function subdirfext ***********/
+ char *subdirfext(char fileres[], char *preop, char *postop)
+ {
+   strcpy(tmpout,preop);
+   strcat(tmpout,fileres);
+   strcat(tmpout,postop);
+   return tmpout;
+ }
+ /*************** function subdirfext3 ***********/
+ char *subdirfext3(char fileres[], char *preop, char *postop)
+ {
+   /* Caution optionfilefiname is hidden */
+   strcpy(tmpout,optionfilefiname);
+   strcat(tmpout,"/");
+   strcat(tmpout,preop);
+   strcat(tmpout,fileres);
+   strcat(tmpout,postop);
+   return tmpout;
+ }
  char *asc_diff_time(long time_sec, char ascdiff[])
  {
    long sec_left, days, hours, minutes;
- Line 1977  void powell(double p[], double **xi, int
+ Line 2059  void powell(double p[], double **xi, int
  double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij)
  {
-   /* Computes the prevalence limit in each live state at age x by left multiplying the unit
+   /* Computes the prevalence limit in each live state at age x and for covariate ij by left multiplying the unit
       matrix by transitions matrix until convergence is reached with precision ftolpl */
    /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1  = Wx-n Px-n ... Px-2 Px-1 I */
    /* Wx is row vector: population in state 1, population in state 2, population dead */
- Line 2000  double **prevalim(double **prlim, int nl
+ Line 2082  double **prevalim(double **prlim, int nl
    int i, ii,j,k;
    double *min, *max, *meandiff, maxmax,sumnew=0.;
    /* double **matprod2(); */ /* test */
-   double **out, cov[NCOVMAX+1], **pmij();
+   double **out, cov[NCOVMAX+1], **pmij(); /* **pmmij is a global variable feeded with oldms etc */
    double **newm;
    double agefin, delaymax=200. ; /* 100 Max number of years to converge */
    int ncvloop=0;
- Line 2009  double **prevalim(double **prlim, int nl
+ Line 2091  double **prevalim(double **prlim, int nl
    max=vector(1,nlstate);
    meandiff=vector(1,nlstate);
+         /* Starting with matrix unity */
    for (ii=1;ii<=nlstate+ndeath;ii++)
      for (j=1;j<=nlstate+ndeath;j++){
        oldm[ii][j]=(ii==j ? 1.0 : 0.0);
- Line 2027  double **prevalim(double **prlim, int nl
+ Line 2110  double **prevalim(double **prlim, int nl
        cov[3]= agefin*agefin;;
      for (k=1; k<=cptcovn;k++) {
        /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+                         /* Here comes the value of the covariate 'ij' */
        cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
        /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */
      }
- Line 2042  double **prevalim(double **prlim, int nl
+ Line 2126  double **prevalim(double **prlim, int nl
      /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
      /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
      /* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */
+                 /* age and covariate values of ij are in 'cov' */
      out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /* Bug Valgrind */
      savm=oldm;
- Line 2088  Earliest age to start was %d-%d=%d, ncvl
+ Line 2173  Earliest age to start was %d-%d=%d, ncvl
    return prlim; /* should not reach here */
  }
+  /**** Back Prevalence limit (stable or period prevalence)  ****************/
+  /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ageminpar, double agemaxpar, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */
+  /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */
+  double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ftolpl, int *ncvyear, int ij)
+ {
+   /* Computes the prevalence limit in each live state at age x and covariate ij by left multiplying the unit
+      matrix by transitions matrix until convergence is reached with precision ftolpl */
+   /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1  = Wx-n Px-n ... Px-2 Px-1 I */
+   /* Wx is row vector: population in state 1, population in state 2, population dead */
+   /* or prevalence in state 1, prevalence in state 2, 0 */
+   /* newm is the matrix after multiplications, its rows are identical at a factor */
+   /* Initial matrix pimij */
+   /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */
+   /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */
+   /*  0,                   0                  , 1} */
+   /*
+    * and after some iteration: */
+   /* {0.45504275246439968, 0.42731458730878791, 0.11764266022681241, */
+   /*  0.45201005341706885, 0.42865420071559901, 0.11933574586733192, */
+   /*  0,                   0                  , 1} */
+   /* And prevalence by suppressing the deaths are close to identical rows in prlim: */
+   /* {0.51571254859325999, 0.4842874514067399, */
+   /*  0.51326036147820708, 0.48673963852179264} */
+   /* If we start from prlim again, prlim tends to a constant matrix */
+   int i, ii,j,k;
+   double *min, *max, *meandiff, maxmax,sumnew=0.;
+   /* double **matprod2(); */ /* test */
+   double **out, cov[NCOVMAX+1], **bmij();
+   double **newm;
+   double         **dnewm, **doldm, **dsavm;  /* for use */
+   double         **oldm, **savm;  /* for use */
+   double agefin, delaymax=200. ; /* 100 Max number of years to converge */
+   int ncvloop=0;
+   min=vector(1,nlstate);
+   max=vector(1,nlstate);
+   meandiff=vector(1,nlstate);
+         dnewm=ddnewms; doldm=ddoldms; dsavm=ddsavms;
+         oldm=oldms; savm=savms;
+         /* Starting with matrix unity */
+         for (ii=1;ii<=nlstate+ndeath;ii++)
+                 for (j=1;j<=nlstate+ndeath;j++){
+       oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+     }
+   cov[1]=1.;
+   /* Even if hstepm = 1, at least one multiplication by the unit matrix */
+   /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */
+   /* for(agefin=age+stepm/YEARM; agefin<=age+delaymax; agefin=agefin+stepm/YEARM){ /\* A changer en age *\/ */
+   for(agefin=age; agefin<AGESUP; agefin=agefin+stepm/YEARM){ /* A changer en age */
+     ncvloop++;
+     newm=savm; /* oldm should be kept from previous iteration or unity at start */
+                 /* newm points to the allocated table savm passed by the function it can be written, savm could be reallocated */
+     /* Covariates have to be included here again */
+     cov[2]=agefin;
+     if(nagesqr==1)
+       cov[3]= agefin*agefin;;
+     for (k=1; k<=cptcovn;k++) {
+       /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+       cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+       /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */
+     }
+     /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+     /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */
+     for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2];
+     for (k=1; k<=cptcovprod;k++) /* Useless */
+       /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+       cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
+     /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
+     /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
+     /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
+     /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
+     /* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */
+                 /* ij should be linked to the correct index of cov */
+                 /* age and covariate values ij are in 'cov', but we need to pass
+                  * ij for the observed prevalence at age and status and covariate
+                  * number:  prevacurrent[(int)agefin][ii][ij]
+                  */
+     /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, ageminpar, agemaxpar, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */
+     /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */
+     out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij)); /* Bug Valgrind */
+     savm=oldm;
+     oldm=newm;
+     for(j=1; j<=nlstate; j++){
+       max[j]=0.;
+       min[j]=1.;
+     }
+     for(j=1; j<=nlstate; j++){
+       for(i=1;i<=nlstate;i++){
+                                 /* bprlim[i][j]= newm[i][j]/(1-sumnew); */
+                                 bprlim[i][j]= newm[i][j];
+                                 max[i]=FMAX(max[i],bprlim[i][j]); /* Max in line */
+                                 min[i]=FMIN(min[i],bprlim[i][j]);
+       }
+     }
+     maxmax=0.;
+     for(i=1; i<=nlstate; i++){
+       meandiff[i]=(max[i]-min[i])/(max[i]+min[i])*2.; /* mean difference for each column */
+       maxmax=FMAX(maxmax,meandiff[i]);
+       /* printf("Back age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, i, meandiff[i],(int)agefin, i, max[i], i, min[i],maxmax); */
+     } /* j loop */
+     *ncvyear= -( (int)age- (int)agefin);
+     /* printf("Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear);*/
+     if(maxmax < ftolpl){
+       printf("OK Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear);
+       free_vector(min,1,nlstate);
+       free_vector(max,1,nlstate);
+       free_vector(meandiff,1,nlstate);
+       return bprlim;
+     }
+   } /* age loop */
+     /* After some age loop it doesn't converge */
+   printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
+ Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
+   /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
+   free_vector(min,1,nlstate);
+   free_vector(max,1,nlstate);
+   free_vector(meandiff,1,nlstate);
+   return bprlim; /* should not reach here */
+ }
  /*************** transition probabilities ***************/
  double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
- Line 2109  double **pmij(double **ps, double *cov,
+ Line 2325  double **pmij(double **ps, double *cov,
    /*double t34;*/
    int i,j, nc, ii, jj;
-     for(i=1; i<= nlstate; i++){
+         for(i=1; i<= nlstate; i++){
-       for(j=1; j<i;j++){
+                 for(j=1; j<i;j++){
-         for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+                         for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
-           /*lnpijopii += param[i][j][nc]*cov[nc];*/
+                                 /*lnpijopii += param[i][j][nc]*cov[nc];*/
-           lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
+                                 lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
- /*       printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+                                 /*       printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
-         }
+                         }
-         ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+                         ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
- /*      printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+                         /*      printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
-       }
+                 }
-       for(j=i+1; j<=nlstate+ndeath;j++){
+                 for(j=i+1; j<=nlstate+ndeath;j++){
-         for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+                         for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
-           /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
+                                 /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
-           lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
+                                 lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
- /*        printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
+                                 /*        printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
-         }
+                         }
-         ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+                         ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
-       }
+                 }
-     }
+         }
-     for(i=1; i<= nlstate; i++){
+         for(i=1; i<= nlstate; i++){
-       s1=0;
+                 s1=0;
-       for(j=1; j<i; j++){
+                 for(j=1; j<i; j++){
-         s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+                         s1+=exp(ps[i][j]); /* In fact sums pij/pii */
-         /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+                         /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
-       }
+                 }
-       for(j=i+1; j<=nlstate+ndeath; j++){
+                 for(j=i+1; j<=nlstate+ndeath; j++){
-         s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+                         s1+=exp(ps[i][j]); /* In fact sums pij/pii */
-         /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+                         /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
-       }
+                 }
-       /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
+                 /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
-       ps[i][i]=1./(s1+1.);
+                 ps[i][i]=1./(s1+1.);
-       /* Computing other pijs */
+                 /* Computing other pijs */
-       for(j=1; j<i; j++)
+                 for(j=1; j<i; j++)
-         ps[i][j]= exp(ps[i][j])*ps[i][i];
+                         ps[i][j]= exp(ps[i][j])*ps[i][i];
-       for(j=i+1; j<=nlstate+ndeath; j++)
+                 for(j=i+1; j<=nlstate+ndeath; j++)
-         ps[i][j]= exp(ps[i][j])*ps[i][i];
+                         ps[i][j]= exp(ps[i][j])*ps[i][i];
-       /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
+                 /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
-     } /* end i */
+         } /* end i */
-     for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
+         for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
-       for(jj=1; jj<= nlstate+ndeath; jj++){
+                 for(jj=1; jj<= nlstate+ndeath; jj++){
-         ps[ii][jj]=0;
+                         ps[ii][jj]=0;
-         ps[ii][ii]=1;
+                         ps[ii][ii]=1;
-       }
+                 }
-     }
+         }
-     /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
+         /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
-     /*   for(jj=1; jj<= nlstate+ndeath; jj++){ */
+         /*   for(jj=1; jj<= nlstate+ndeath; jj++){ */
-     /*  printf(" pmij  ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
+         /*      printf(" pmij  ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
-     /*   } */
+         /*   } */
-     /*   printf("\n "); */
+         /*   printf("\n "); */
-     /* } */
+         /* } */
-     /* printf("\n ");printf("%lf ",cov[2]);*/
+         /* printf("\n ");printf("%lf ",cov[2]);*/
-     /*
+         /*
-       for(i=1; i<= npar; i++) printf("%f ",x[i]);
+                 for(i=1; i<= npar; i++) printf("%f ",x[i]);
-       goto end;*/
+                 goto end;*/
-     return ps;
+         return ps;
+ }
+ /*************** backward transition probabilities ***************/
+  /* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate,  double ***prevacurrent, double ageminpar, double agemaxpar, double ***dnewm, double **doldm, double **dsavm, int ij ) */
+ /* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate,  double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */
+  double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate,  double ***prevacurrent, int ij )
+ {
+         /* Computes the backward probability at age agefin and covariate ij
+          * and returns in **ps as well as **bmij.
+          */
+   int i, ii, j,k;
+         double **out, **pmij();
+         double sumnew=0.;
+   double agefin;
+         double **dnewm, **dsavm, **doldm;
+         double **bbmij;
+   doldm=ddoldms; /* global pointers */
+         dnewm=ddnewms;
+         dsavm=ddsavms;
+         agefin=cov[2];
+         /* bmij *//* age is cov[2], ij is included in cov, but we need for
+                  the observed prevalence (with this covariate ij) */
+         dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate);
+         /* We do have the matrix Px in savm  and we need pij */
+         for (j=1;j<=nlstate+ndeath;j++){
+                 sumnew=0.; /* w1 p11 + w2 p21 only on live states */
+                 for (ii=1;ii<=nlstate;ii++){
+                         sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij];
+                 } /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */
+                 for (ii=1;ii<=nlstate+ndeath;ii++){
+                         if(sumnew >= 1.e-10){
+                                 /* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */
+                                 /*      doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
+                                 /* }else if(agefin >= agemaxpar+stepm/YEARM){ */
+                                 /*      doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
+                                 /* }else */
+                                         doldm[ii][j]=(ii==j ? 1./sumnew : 0.0);
+                         }else{
+                                 printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin);
+                         }
+                 } /*End ii */
+         } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */
+                 /* left Product of this diag matrix by dsavm=Px (newm=dsavm*doldm) */
+         bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */
+         /* dsavm=doldm; /\* dsavm is now diag [1/(w_1p1i+w_2 p2i)] but can be overwritten*\/ */
+         /* doldm=dnewm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */
+         /* dnewm=dsavm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */
+         /* left Product of this matrix by diag matrix of prevalences (savm) */
+         for (j=1;j<=nlstate+ndeath;j++){
+                 for (ii=1;ii<=nlstate+ndeath;ii++){
+                         dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0);
+                 }
+         } /* End j, At the end oldm is diag[1/(w_1p1i+w_2 p2i)] */
+         ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */
+         /* newm or out is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */
+         /* end bmij */
+         return ps;
+ }
+ /*************** transition probabilities ***************/
+ double **bpmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
+ {
+   /* According to parameters values stored in x and the covariate's values stored in cov,
+      computes the probability to be observed in state j being in state i by appying the
+      model to the ncovmodel covariates (including constant and age).
+      lnpijopii=ln(pij/pii)= aij+bij*age+cij*v1+dij*v2+... = sum_nc=1^ncovmodel xij(nc)*cov[nc]
+      and, according on how parameters are entered, the position of the coefficient xij(nc) of the
+      ncth covariate in the global vector x is given by the formula:
+      j<i nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel
+      j>=i nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel
+      Computes ln(pij/pii) (lnpijopii), deduces pij/pii by exponentiation,
+      sums on j different of i to get 1-pii/pii, deduces pii, and then all pij.
+      Outputs ps[i][j] the probability to be observed in j being in j according to
+      the values of the covariates cov[nc] and corresponding parameter values x[nc+shiftij]
+   */
+   double s1, lnpijopii;
+   /*double t34;*/
+   int i,j, nc, ii, jj;
+         for(i=1; i<= nlstate; i++){
+                 for(j=1; j<i;j++){
+                         for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+                                 /*lnpijopii += param[i][j][nc]*cov[nc];*/
+                                 lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
+                                 /*       printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+                         }
+                         ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+                         /*      printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+                 }
+                 for(j=i+1; j<=nlstate+ndeath;j++){
+                         for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+                                 /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
+                                 lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
+                                 /*        printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
+                         }
+                         ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+                 }
+         }
+         for(i=1; i<= nlstate; i++){
+                 s1=0;
+                 for(j=1; j<i; j++){
+                         s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+                         /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+                 }
+                 for(j=i+1; j<=nlstate+ndeath; j++){
+                         s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+                         /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+                 }
+                 /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
+                 ps[i][i]=1./(s1+1.);
+                 /* Computing other pijs */
+                 for(j=1; j<i; j++)
+                         ps[i][j]= exp(ps[i][j])*ps[i][i];
+                 for(j=i+1; j<=nlstate+ndeath; j++)
+                         ps[i][j]= exp(ps[i][j])*ps[i][i];
+                 /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
+         } /* end i */
+         for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
+                 for(jj=1; jj<= nlstate+ndeath; jj++){
+                         ps[ii][jj]=0;
+                         ps[ii][ii]=1;
+                 }
+         }
+         /* Added for backcast */ /* Transposed matrix too */
+         for(jj=1; jj<= nlstate+ndeath; jj++){
+                 s1=0.;
+                 for(ii=1; ii<= nlstate+ndeath; ii++){
+                         s1+=ps[ii][jj];
+                 }
+                 for(ii=1; ii<= nlstate; ii++){
+                         ps[ii][jj]=ps[ii][jj]/s1;
+                 }
+         }
+         /* Transposition */
+         for(jj=1; jj<= nlstate+ndeath; jj++){
+                 for(ii=jj; ii<= nlstate+ndeath; ii++){
+                         s1=ps[ii][jj];
+                         ps[ii][jj]=ps[jj][ii];
+                         ps[jj][ii]=s1;
+                 }
+         }
+         /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
+         /*   for(jj=1; jj<= nlstate+ndeath; jj++){ */
+         /*      printf(" pmij  ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
+         /*   } */
+         /*   printf("\n "); */
+         /* } */
+         /* printf("\n ");printf("%lf ",cov[2]);*/
+         /*
+                 for(i=1; i<= npar; i++) printf("%f ",x[i]);
+                 goto end;*/
+         return ps;
  }
  /**************** Product of 2 matrices ******************/
  double **matprod2(double **out, double **in,int nrl, int nrh, int ncl, int nch, int ncolol, int ncoloh, double **b)
- Line 2194  double **matprod2(double **out, double *
+ Line 2570  double **matprod2(double **out, double *
  double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij )
  {
-   /* Computes the transition matrix starting at age 'age' over
+   /* Computes the transition matrix starting at age 'age' and combination of covariate values corresponding to ij over
       'nhstepm*hstepm*stepm' months (i.e. until
       age (in years)  age+nhstepm*hstepm*stepm/12) by multiplying
       nhstepm*hstepm matrices.
- Line 2210  double ***hpxij(double ***po, int nhstep
+ Line 2586  double ***hpxij(double ***po, int nhstep
    double **out, cov[NCOVMAX+1];
    double **newm;
    double agexact;
+   double agebegin, ageend;
    /* Hstepm could be zero and should return the unit matrix */
    for (i=1;i<=nlstate+ndeath;i++)
- Line 2223  double ***hpxij(double ***po, int nhstep
+ Line 2600  double ***hpxij(double ***po, int nhstep
        newm=savm;
        /* Covariates have to be included here again */
        cov[1]=1.;
-       agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
+       agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */
        cov[2]=agexact;
        if(nagesqr==1)
-         cov[3]= agexact*agexact;
+                                 cov[3]= agexact*agexact;
        for (k=1; k<=cptcovn;k++)
-         cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+                                 cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
-         /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+                         /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
        for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
-         /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+                                 /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-         cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+                                 cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
-         /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
+                         /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
        for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
-         cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
+                                 cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
-         /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+                         /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
        /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
        /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
+                         /* right multiplication of oldm by the current matrix */
        out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,
                     pmij(pmmij,cov,ncovmodel,x,nlstate));
+       /* if((int)age == 70){ */
+       /*        printf(" Forward hpxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */
+       /*        for(i=1; i<=nlstate+ndeath; i++) { */
+       /*          printf("%d pmmij ",i); */
+       /*          for(j=1;j<=nlstate+ndeath;j++) { */
+       /*            printf("%f ",pmmij[i][j]); */
+       /*          } */
+       /*          printf(" oldm "); */
+       /*          for(j=1;j<=nlstate+ndeath;j++) { */
+       /*            printf("%f ",oldm[i][j]); */
+       /*          } */
+       /*          printf("\n"); */
+       /*        } */
+       /* } */
        savm=oldm;
        oldm=newm;
      }
      for(i=1; i<=nlstate+ndeath; i++)
        for(j=1;j<=nlstate+ndeath;j++) {
-         po[i][j][h]=newm[i][j];
+                                 po[i][j][h]=newm[i][j];
-         /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
+                                 /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
        }
      /*printf("h=%d ",h);*/
    } /* end h */
- /*     printf("\n H=%d \n",h); */
+         /*     printf("\n H=%d \n",h); */
    return po;
  }
+ /************* Higher Back Matrix Product ***************/
+ /* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */
+  double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij )
+ {
+   /* Computes the transition matrix starting at age 'age' over
+      'nhstepm*hstepm*stepm' months (i.e. until
+      age (in years)  age+nhstepm*hstepm*stepm/12) by multiplying
+      nhstepm*hstepm matrices.
+      Output is stored in matrix po[i][j][h] for h every 'hstepm' step
+      (typically every 2 years instead of every month which is too big
+      for the memory).
+      Model is determined by parameters x and covariates have to be
+      included manually here.
+      */
+   int i, j, d, h, k;
+   double **out, cov[NCOVMAX+1];
+   double **newm;
+   double agexact;
+   double agebegin, ageend;
+         double **oldm, **savm;
+         oldm=oldms;savm=savms;
+   /* Hstepm could be zero and should return the unit matrix */
+   for (i=1;i<=nlstate+ndeath;i++)
+     for (j=1;j<=nlstate+ndeath;j++){
+       oldm[i][j]=(i==j ? 1.0 : 0.0);
+       po[i][j][0]=(i==j ? 1.0 : 0.0);
+     }
+   /* Even if hstepm = 1, at least one multiplication by the unit matrix */
+   for(h=1; h <=nhstepm; h++){
+     for(d=1; d <=hstepm; d++){
+       newm=savm;
+       /* Covariates have to be included here again */
+       cov[1]=1.;
+       agexact=age-((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */
+       /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */
+       cov[2]=agexact;
+       if(nagesqr==1)
+                                 cov[3]= agexact*agexact;
+       for (k=1; k<=cptcovn;k++)
+                                 cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+                         /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+       for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
+                                 /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+                                 cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+                         /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
+       for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
+                                 cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
+                         /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+       /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
+       /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
+       /* Careful transposed matrix */
+                         /* age is in cov[2] */
+       /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */
+                         /*                                               1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */
+       out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\
+,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
+       /* if((int)age == 70){ */
+       /*        printf(" Backward hbxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */
+       /*        for(i=1; i<=nlstate+ndeath; i++) { */
+       /*          printf("%d pmmij ",i); */
+       /*          for(j=1;j<=nlstate+ndeath;j++) { */
+       /*            printf("%f ",pmmij[i][j]); */
+       /*          } */
+       /*          printf(" oldm "); */
+       /*          for(j=1;j<=nlstate+ndeath;j++) { */
+       /*            printf("%f ",oldm[i][j]); */
+       /*          } */
+       /*          printf("\n"); */
+       /*        } */
+       /* } */
+       savm=oldm;
+       oldm=newm;
+     }
+     for(i=1; i<=nlstate+ndeath; i++)
+       for(j=1;j<=nlstate+ndeath;j++) {
+                                 po[i][j][h]=newm[i][j];
+                                 /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
+       }
+     /*printf("h=%d ",h);*/
+   } /* end h */
+         /*     printf("\n H=%d \n",h); */
+   return po;
+ }
  #ifdef NLOPT
    double  myfunc(unsigned n, const double *p1, double *grad, void *pd){
    double fret;
- Line 2400  double func( double *x)
+ Line 2882  double func( double *x)
  /*        else */
  /*          lli=log(out[s1][s2] - savm[s1][s2]); */
  /* #endif */
-             lli=log(out[s1][s2] - savm[s1][s2]);
+           lli=log(out[s1][s2] - savm[s1][s2]);
-         } else if  (s2==-2) {
+         } else if  ( s2==-1 ) { /* alive */
            for (j=1,survp=0. ; j<=nlstate; j++)
              survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
            /*survp += out[s1][j]; */
            lli= log(survp);
          }
          else if  (s2==-4) {
            for (j=3,survp=0. ; j<=nlstate; j++)
              survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
            lli= log(survp);
          }
          else if  (s2==-5) {
            for (j=1,survp=0. ; j<=2; j++)
              survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
            lli= log(survp);
          }
          else{
            lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
            /*  lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */
- Line 2532  double func( double *x)
+ Line 3011  double func( double *x)
          s2=s[mw[mi+1][i]][i];
          if( s2 > nlstate){
            lli=log(out[s1][s2] - savm[s1][s2]);
+         } else if  ( s2==-1 ) { /* alive */
+           for (j=1,survp=0. ; j<=nlstate; j++)
+             survp += out[s1][j];
+           lli= log(survp);
          }else{
            lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
          }
- Line 2594  double funcone( double *x)
+ Line 3077  double funcone( double *x)
    int s1, s2;
    double bbh, survp;
    double agexact;
+   double agebegin, ageend;
    /*extern weight */
    /* We are differentiating ll according to initial status */
    /*  for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
- Line 2612  double funcone( double *x)
+ Line 3096  double funcone( double *x)
            oldm[ii][j]=(ii==j ? 1.0 : 0.0);
            savm[ii][j]=(ii==j ? 1.0 : 0.0);
          }
-       for(d=0; d<dh[mi][i]; d++){
+       agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */
+       ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */
+       for(d=0; d<dh[mi][i]; d++){  /* Delay between two effective waves */
+         /*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
+           and mw[mi+1][i]. dh depends on stepm.*/
          newm=savm;
          agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
          cov[2]=agexact;
- Line 2633  double funcone( double *x)
+ Line 3122  double funcone( double *x)
        s1=s[mw[mi][i]][i];
        s2=s[mw[mi+1][i]][i];
+       /* if(s2==-1){ */
+       /*        printf(" s1=%d, s2=%d i=%d \n", s1, s2, i); */
+       /*        /\* exit(1); *\/ */
+       /* } */
        bbh=(double)bh[mi][i]/(double)stepm;
        /* bias is positive if real duration
         * is higher than the multiple of stepm and negative otherwise.
         */
        if( s2 > nlstate && (mle <5) ){  /* Jackson */
          lli=log(out[s1][s2] - savm[s1][s2]);
-       } else if  (s2==-2) {
+       } else if  ( s2==-1 ) { /* alive */
          for (j=1,survp=0. ; j<=nlstate; j++)
            survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
          lli= log(survp);
- Line 2660  double funcone( double *x)
+ Line 3153  double funcone( double *x)
        ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
        /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
        if(globpr){
-         fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\
+         fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\
   %11.6f %11.6f %11.6f ", \
-                 num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
+                 num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
 *weight[i]*lli,out[s1][s2],savm[s1][s2]);
          for(k=1,llt=0.,l=0.; k<=nlstate; k++){
            llt +=ll[k]*gipmx/gsw;
- Line 2700  void likelione(FILE *ficres,double p[],
+ Line 3193  void likelione(FILE *ficres,double p[],
        printf("Problem with resultfile: %s\n", fileresilk);
        fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);
      }
-     fprintf(ficresilk, "#individual(line's_record) count age s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
+     fprintf(ficresilk, "#individual(line's_record) count ageb ageend s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
-     fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav ");
+     fprintf(ficresilk, "#num_i ageb agend i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav ");
      /*  i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */
      for(k=1; k<=nlstate; k++)
        fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);
- Line 2719  void likelione(FILE *ficres,double p[],
+ Line 3212  void likelione(FILE *ficres,double p[],
      for (k=1; k<= nlstate ; k++) {
-       fprintf(fichtm,"<br>- Probability p%dj by origin %d and destination j <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. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
  <img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k);
      }
      fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \
- Line 3018  double hessij( double x[], double **hess
+ Line 3511  double hessij( double x[], double **hess
        kmax=kmax+10;
      }
      if(kmax >=10 || firstime ==1){
-       printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol);
+       printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol);
-       fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; 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 may 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);
        fprintf(ficlog,"%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 3175  void pstamp(FILE *fichier)
+ Line 3668  void pstamp(FILE *fichier)
  }
  /************ Frequencies ********************/
- void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[])
+ void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
+                   int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[],\
+                   int firstpass,  int lastpass, int stepm, int weightopt, char model[])
  {  /* Some frequencies */
    int i, m, jk, j1, bool, z1,j;
+   int mi; /* Effective wave */
    int first;
    double ***freq; /* Frequencies */
    double *pp, **prop;
    double pos,posprop, k2, dateintsum=0,k2cpt=0;
-   char fileresp[FILENAMELENGTH];
+   char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH];
+   double agebegin, ageend;
    pp=vector(1,nlstate);
-   prop=matrix(1,nlstate,iagemin,iagemax+3);
+   prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
+   /* prop=matrix(1,nlstate,iagemin,iagemax+3); */
    strcpy(fileresp,"P_");
    strcat(fileresp,fileresu);
+   /*strcat(fileresphtm,fileresu);*/
    if((ficresp=fopen(fileresp,"w"))==NULL) {
      printf("Problem with prevalence resultfile: %s\n", fileresp);
      fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
      exit(0);
    }
-   freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3);
+   strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));
+   if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {
+     printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
+     fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
+     fflush(ficlog);
+     exit(70);
+   }
+   else{
+     fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
+ <hr size=\"2\" color=\"#EC5E5E\"> \n\
+ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
+           fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
+   }
+     fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\n",fileresphtm, fileresphtm);
+   strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
+   if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
+     printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
+     fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
+     fflush(ficlog);
+     exit(70);
+   }
+   else{
+     fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
+ <hr size=\"2\" color=\"#EC5E5E\"> \n\
+ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
+           fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
+   }
+   fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
+   freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
    j1=0;
    j=cptcoveff;
- Line 3202  void  freqsummary(char fileres[], int ia
+ Line 3732  void  freqsummary(char fileres[], int ia
    first=1;
-   /* for(k1=1; k1<=j ; k1++){ */  /* Loop on covariates */
+   for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ /* Loop on covariates combination */
-   /*  for(i1=1; i1<=ncodemax[k1];i1++){ */ /* Now it is 2 */
-   /*    j1++; */
-   for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){
        /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
          scanf("%d", i);*/
        for (i=-5; i<=nlstate+ndeath; i++)
- Line 3219  void  freqsummary(char fileres[], int ia
+ Line 3746  void  freqsummary(char fileres[], int ia
        dateintsum=0;
        k2cpt=0;
-       for (i=1; i<=imx; i++) {
+       for (i=1; i<=imx; i++) { /* For each individual i */
          bool=1;
-         if  (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
+         if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
            for (z1=1; z1<=cptcoveff; z1++)
              if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
                  /* Tests if the value of each of the covariates of i is equal to filter j1 */
- Line 3231  void  freqsummary(char fileres[], int ia
+ Line 3758  void  freqsummary(char fileres[], int ia
                  j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
                /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
              }
-         }
+         } /* cptcovn > 0 */
          if (bool==1){
-           for(m=firstpass; m<=lastpass; m++){
+           /* for(m=firstpass; m<=lastpass; m++){ */
-             k2=anint[m][i]+(mint[m][i]/12.);
+           for(mi=1; mi<wav[i];mi++){
-             /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
+             m=mw[mi][i];
-               if(agev[m][i]==0) agev[m][i]=iagemax+1;
+             /* dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective (mi) waves m=mw[mi][i]
-               if(agev[m][i]==1) agev[m][i]=iagemax+2;
+                and mw[mi+1][i]. dh depends on stepm. */
-               if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i];
+             agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
+             ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /* Age at end of wave and transition */
+             if(m >=firstpass && m <=lastpass){
+               k2=anint[m][i]+(mint[m][i]/12.);
+               /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
+               if(agev[m][i]==0) agev[m][i]=iagemax+1;  /* All ages equal to 0 are in iagemax+1 */
+               if(agev[m][i]==1) agev[m][i]=iagemax+2;  /* All ages equal to 1 are in iagemax+2 */
+               if (s[m][i]>0 && s[m][i]<=nlstate)  /* If status at wave m is known and a live state */
+                 prop[s[m][i]][(int)agev[m][i]] += weight[i];  /* At age of beginning of transition, where status is known */
                if (m<lastpass) {
-                 freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
+                 /* if(s[m][i]==4 && s[m+1][i]==4) */
-                 freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i];
+                 /*   printf(" num=%ld m=%d, i=%d s1=%d s2=%d agev at m=%d\n", num[i], m, i,s[m][i],s[m+1][i], (int)agev[m][i]); */
+                 if(s[m][i]==-1)
+                   printf(" num=%ld m=%d, i=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[i], m, i,s[m][i],s[m+1][i], (int)agev[m][i],agebegin, ageend, (int)((agebegin+ageend)/2.));
+                 freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; /* At age of beginning of transition, where status is known */
+                 /* freq[s[m][i]][s[m+1][i]][(int)((agebegin+ageend)/2.)] += weight[i]; */
+                 freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
                }
+             }
-               if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3))) {
+             if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3)) && (anint[m][i]!=9999) && (mint[m][i]!=99)) {
-                 dateintsum=dateintsum+k2;
+               dateintsum=dateintsum+k2;
-                 k2cpt++;
+               k2cpt++;
-               }
+               /* printf("i=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",i, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
-               /*}*/
+             }
-           }
+             /*}*/
-         }
+           } /* end m */
-       } /* end i */
+         } /* end bool */
+       } /* end i = 1 to imx */
        /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
        pstamp(ficresp);
        if  (cptcovn>0) {
          fprintf(ficresp, "\n#********** Variable ");
-         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+         fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");
-         fprintf(ficresp, "**********\n#");
+         fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");
+         for (z1=1; z1<=cptcoveff; z1++){
+           fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+           fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+           fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+         }
+           fprintf(ficresp, "**********\n#");
+         fprintf(ficresphtm, "**********</h3>\n");
+         fprintf(ficresphtmfr, "**********</h3>\n");
          fprintf(ficlog, "\n#********** Variable ");
          for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-         fprintf(ficlog, "**********\n#");
+         fprintf(ficlog, "**********\n");
        }
-       for(i=1; i<=nlstate;i++)
+       fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
+       for(i=1; i<=nlstate;i++) {
          fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
+         fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);
+       }
        fprintf(ficresp, "\n");
+       fprintf(ficresphtm, "\n");
+       /* Header of frequency table by age */
+       fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");
+       fprintf(ficresphtmfr,"<th>Age</th> ");
+       for(jk=-1; jk <=nlstate+ndeath; jk++){
+         for(m=-1; m <=nlstate+ndeath; m++){
+           if(jk!=0 && m!=0)
+             fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);
+         }
+       }
+       fprintf(ficresphtmfr, "\n");
+       /* For each age */
        for(i=iagemin; i <= iagemax+3; i++){
-         if(i==iagemax+3){
+         fprintf(ficresphtm,"<tr>");
+         if(i==iagemax+1){
+           fprintf(ficlog,"1");
+           fprintf(ficresphtmfr,"<tr><th>0</th> ");
+         }else if(i==iagemax+2){
+           fprintf(ficlog,"0");
+           fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
+         }else if(i==iagemax+3){
            fprintf(ficlog,"Total");
+           fprintf(ficresphtmfr,"<tr><th>Total</th> ");
          }else{
            if(first==1){
              first=0;
              printf("See log file for details...\n");
            }
+           fprintf(ficresphtmfr,"<tr><th>%d</th> ",i);
            fprintf(ficlog,"Age %d", i);
          }
          for(jk=1; jk <=nlstate ; jk++){
- Line 3318  void  freqsummary(char fileres[], int ia
+ Line 3892  void  freqsummary(char fileres[], int ia
            if( i <= iagemax){
              if(pos>=1.e-5){
                fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop);
+               fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",i,prop[jk][i]/posprop, prop[jk][i],posprop);
                /*probs[i][jk][j1]= pp[jk]/pos;*/
                /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/
              }
-             else
+             else{
                fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop);
+               fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",i, prop[jk][i],posprop);
+             }
            }
          }
-         for(jk=-1; jk <=nlstate+ndeath; jk++)
+         for(jk=-1; jk <=nlstate+ndeath; jk++){
-           for(m=-1; m <=nlstate+ndeath; m++)
+           for(m=-1; m <=nlstate+ndeath; m++){
-             if(freq[jk][m][i] !=0 ) {
+             if(freq[jk][m][i] !=0 ) { /* minimizing output */
-             if(first==1)
+               if(first==1){
-               printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
+                 printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
+               }
                fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]);
              }
-         if(i <= iagemax)
+             if(jk!=0 && m!=0)
+               fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][i]);
+           }
+         }
+         fprintf(ficresphtmfr,"</tr>\n ");
+         if(i <= iagemax){
            fprintf(ficresp,"\n");
+           fprintf(ficresphtm,"</tr>\n");
+         }
          if(first==1)
            printf("Others in log...\n");
          fprintf(ficlog,"\n");
-       }
+       } /* end loop i */
+       fprintf(ficresphtm,"</table>\n");
+       fprintf(ficresphtmfr,"</table>\n");
        /*}*/
-   }
+   } /* end j1 */
    dateintmean=dateintsum/k2cpt;
    fclose(ficresp);
-   free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3);
+   fclose(ficresphtm);
+   fclose(ficresphtmfr);
+   free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+3+AGEMARGE);
    free_vector(pp,1,nlstate);
-   free_matrix(prop,1,nlstate,iagemin, iagemax+3);
+   free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE);
    /* End of Freq */
  }
- Line 3359  void prevalence(double ***probs, double
+ Line 3948  void prevalence(double ***probs, double
    */
    int i, m, jk, j1, bool, z1,j;
+   int mi; /* Effective wave */
+   int iage;
+   double agebegin, ageend;
    double **prop;
    double posprop;
- Line 3369  void prevalence(double ***probs, double
+ Line 3961  void prevalence(double ***probs, double
    iagemin= (int) agemin;
    iagemax= (int) agemax;
    /*pp=vector(1,nlstate);*/
-   prop=matrix(1,nlstate,iagemin,iagemax+3);
+   prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
    /*  freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/
    j1=0;
- Line 3378  void prevalence(double ***probs, double
+ Line 3970  void prevalence(double ***probs, double
    first=1;
    for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){
-     /*for(i1=1; i1<=ncodemax[k1];i1++){
+     for (i=1; i<=nlstate; i++)
-       j1++;*/
+       for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++)
+                                 prop[i][iage]=0.0;
+     for (i=1; i<=imx; i++) { /* Each individual */
+       bool=1;
+       if  (cptcovn>0) {  /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
+                                 for (z1=1; z1<=cptcoveff; z1++)
+                                         if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)])
+                                                 bool=0;
+       }
+       if (bool==1) {
+                                 /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */
+                                 for(mi=1; mi<wav[i];mi++){
+                                         m=mw[mi][i];
+                                         agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
+                                         /* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */
+                                         if(m >=firstpass && m <=lastpass){
+                                                 y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */
+                                                 if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
+                                                         if(agev[m][i]==0) agev[m][i]=iagemax+1;
+                                                         if(agev[m][i]==1) agev[m][i]=iagemax+2;
+                                                         if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+3+AGEMARGE){
+                                                                 printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d  m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m);
+                                                                 exit(1);
+                                                         }
+                                                         if (s[m][i]>0 && s[m][i]<=nlstate) {
+                                                                 /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/
+                                                                 prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */
+                                                                 prop[s[m][i]][iagemax+3] += weight[i];
+                                                         } /* end valid statuses */
+                                                 } /* end selection of dates */
+                                         } /* end selection of waves */
+                                 } /* end effective waves */
+       } /* end bool */
+     }
+     for(i=iagemin; i <= iagemax+3; i++){
+       for(jk=1,posprop=0; jk <=nlstate ; jk++) {
+                                 posprop += prop[jk][i];
+       }
-       for (i=1; i<=nlstate; i++)
+       for(jk=1; jk <=nlstate ; jk++){
-         for(m=iagemin; m <= iagemax+3; m++)
+                                 if( i <=  iagemax){
-           prop[i][m]=0.0;
+                                         if(posprop>=1.e-5){
+                                                 probs[i][jk][j1]= prop[jk][i]/posprop;
-       for (i=1; i<=imx; i++) { /* Each individual */
+                                         } else{
-         bool=1;
+                                                 if(first==1){
-         if  (cptcovn>0) {
+                                                         first=0;
-           for (z1=1; z1<=cptcoveff; z1++)
+                                                         printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]);
-             if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)])
+                                                 }
-               bool=0;
+                                         }
-         }
+                                 }
-         if (bool==1) {
+       }/* end jk */
-           for(m=firstpass; m<=lastpass; m++){/* Other selection (we can limit to certain interviews*/
+     }/* end i */
-             y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */
-             if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
-               if(agev[m][i]==0) agev[m][i]=iagemax+1;
-               if(agev[m][i]==1) agev[m][i]=iagemax+2;
-               if((int)agev[m][i] <iagemin || (int)agev[m][i] >iagemax+3) printf("Error on individual =%d agev[m][i]=%f m=%d\n",i, agev[m][i],m);
-               if (s[m][i]>0 && s[m][i]<=nlstate) {
-                 /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/
-                 prop[s[m][i]][(int)agev[m][i]] += weight[i];
-                 prop[s[m][i]][iagemax+3] += weight[i];
-               }
-             }
-           } /* end selection of waves */
-         }
-       }
-       for(i=iagemin; i <= iagemax+3; i++){
-         for(jk=1,posprop=0; jk <=nlstate ; jk++) {
-           posprop += prop[jk][i];
-         }
-         for(jk=1; jk <=nlstate ; jk++){
-           if( i <=  iagemax){
-             if(posprop>=1.e-5){
-               probs[i][jk][j1]= prop[jk][i]/posprop;
-             } else{
-               if(first==1){
-                 first=0;
-                 printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]);
-               }
-             }
-           }
-         }/* end jk */
-       }/* end i */
      /*} *//* end i1 */
    } /* end j1 */
    /*  free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/
    /*free_vector(pp,1,nlstate);*/
-   free_matrix(prop,1,nlstate, iagemin,iagemax+3);
+   free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE);
  }  /* End of prevalence */
  /************* Waves Concatenation ***************/
- Line 3448  void  concatwav(int wav[], int **dh, int
+ Line 4046  void  concatwav(int wav[], int **dh, int
    int i, mi, m;
    /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
       double sum=0., jmean=0.;*/
-   int first;
+   int first, firstwo, firsthree, firstfour;
    int j, k=0,jk, ju, jl;
    double sum=0.;
    first=0;
+   firstwo=0;
+   firsthree=0;
+   firstfour=0;
    jmin=100000;
    jmax=-1;
    jmean=0.;
-   for(i=1; i<=imx; i++){
+   for(i=1; i<=imx; i++){  /* For simple cases and if state is death */
      mi=0;
      m=firstpass;
-     while(s[m][i] <= nlstate){
+     while(s[m][i] <= nlstate){  /* a live state */
-       if(s[m][i]>=1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5)
+       if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */
          mw[++mi][i]=m;
-       if(m >=lastpass)
+       }
+       if(m >=lastpass){
+         if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){
+           if(firsthree == 0){
+             printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
+             firsthree=1;
+           }
+           fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
+           mw[++mi][i]=m;
+         }
+         if(s[m][i]==-2){ /* Vital status is really unknown */
+           nbwarn++;
+           if((int)anint[m][i] == 9999){  /*  Has the vital status really been verified? */
+             printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
+             fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
+           }
+           break;
+         }
          break;
+       }
        else
          m++;
      }/* end while */
-     if (s[m][i] > nlstate){
+     /* After last pass */
+     if (s[m][i] > nlstate){  /* In a death state */
        mi++;     /* Death is another wave */
        /* if(mi==0)  never been interviewed correctly before death */
           /* Only death is a correct wave */
        mw[mi][i]=m;
+     }else if ((int) andc[i] != 9999) { /* Status is either death or negative. A death occured after lastpass, we can't take it into account because of potential bias */
+       /* m++; */
+       /* mi++; */
+       /* s[m][i]=nlstate+1;  /\* We are setting the status to the last of non live state *\/ */
+       /* mw[mi][i]=m; */
+       nberr++;
+       if ((int)anint[m][i]!= 9999) { /* date of last interview is known */
+         if(firstwo==0){
+           printf("Error! Death for individual %ld line=%d  occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+           firstwo=1;
+         }
+         fprintf(ficlog,"Error! Death for individual %ld line=%d  occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+       }else{ /* end date of interview is known */
+         /* death is known but not confirmed by death status at any wave */
+         if(firstfour==0){
+           printf("Error! Death for individual %ld line=%d  occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+           firstfour=1;
+         }
+         fprintf(ficlog,"Error! Death for individual %ld line=%d  occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+       }
      }
      wav[i]=mi;
      if(mi==0){
        nbwarn++;
- Line 3485  void  concatwav(int wav[], int **dh, int
+ Line 4125  void  concatwav(int wav[], int **dh, int
        }
      } /* end mi==0 */
    } /* End individuals */
+   /* wav and mw are no more changed */
    for(i=1; i<=imx; i++){
      for(mi=1; mi<wav[i];mi++){
        if (stepm <=0)
- Line 3915  void cvevsij(double ***eij, double x[],
+ Line 4557  void cvevsij(double ***eij, double x[],
      /* Typically if 20 years nstepm = 20*12/6=40 stepm */
      /* if (stepm >= YEARM) hstepm=1;*/
      nhstepma = nstepma/hstepm;/* Expressed in hstepm, typically nhstepma=40/4=10 */
      /* If stepm=6 months */
      /* Computed by stepm unit matrices, product of hstepma matrices, stored
         in an array of nhstepma length: nhstepma=10, hstepm=4, stepm=6 months */
      hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
      /* Computing  Variances of health expectancies */
      /* Gradient is computed with plus gp and minus gm. Code is duplicated in order to
         decrease memory allocation */
      for(theta=1; theta <=npar; theta++){
        for(i=1; i<=npar; i++){
-         xp[i] = x[i] + (i==theta ?delti[theta]:0);
+                                 xp[i] = x[i] + (i==theta ?delti[theta]:0);
-         xm[i] = x[i] - (i==theta ?delti[theta]:0);
+                                 xm[i] = x[i] - (i==theta ?delti[theta]:0);
        }
        hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij);
        hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij);
        for(j=1; j<= nlstate; j++){
-         for(i=1; i<=nlstate; i++){
+                                 for(i=1; i<=nlstate; i++){
-           for(h=0; h<=nhstepm-1; h++){
+                                         for(h=0; h<=nhstepm-1; h++){
-             gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.;
+                                                 gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.;
-             gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.;
+                                                 gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.;
-           }
+                                         }
-         }
+                                 }
        }
        for(ij=1; ij<= nlstate*nlstate; ij++)
-         for(h=0; h<=nhstepm-1; h++){
+                                 for(h=0; h<=nhstepm-1; h++){
-           gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta];
+                                         gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta];
-         }
+                                 }
      }/* End theta */
      for(h=0; h<=nhstepm-1; h++)
        for(j=1; j<=nlstate*nlstate;j++)
-         for(theta=1; theta <=npar; theta++)
+                                 for(theta=1; theta <=npar; theta++)
-           trgradg[h][j][theta]=gradg[h][theta][j];
+                                         trgradg[h][j][theta]=gradg[h][theta][j];
-      for(ij=1;ij<=nlstate*nlstate;ij++)
+                 for(ij=1;ij<=nlstate*nlstate;ij++)
        for(ji=1;ji<=nlstate*nlstate;ji++)
-         varhe[ij][ji][(int)age] =0.;
+                                 varhe[ij][ji][(int)age] =0.;
-      printf("%d|",(int)age);fflush(stdout);
+                 printf("%d|",(int)age);fflush(stdout);
-      fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
+                 fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
-      for(h=0;h<=nhstepm-1;h++){
+                 for(h=0;h<=nhstepm-1;h++){
        for(k=0;k<=nhstepm-1;k++){
-         matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov);
+                                 matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov);
-         matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]);
+                                 matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]);
-         for(ij=1;ij<=nlstate*nlstate;ij++)
+                                 for(ij=1;ij<=nlstate*nlstate;ij++)
-           for(ji=1;ji<=nlstate*nlstate;ji++)
+                                         for(ji=1;ji<=nlstate*nlstate;ji++)
-             varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;
+                                                 varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;
        }
      }
      /* Computing expectancies */
      hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij);
      for(i=1; i<=nlstate;i++)
        for(j=1; j<=nlstate;j++)
-         for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
+                                 for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
-           eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[i][j][h+1])/2.0*hf;
+                                         eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[i][j][h+1])/2.0*hf;
-           /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
+                                         /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
-         }
+                                 }
      fprintf(ficresstdeij,"%3.0f",age );
      for(i=1; i<=nlstate;i++){
        eip=0.;
        vip=0.;
        for(j=1; j<=nlstate;j++){
-         eip += eij[i][j][(int)age];
+                                 eip += eij[i][j][(int)age];
-         for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */
+                                 for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */
-           vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age];
+                                         vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age];
-         fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) );
+                                 fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) );
        }
        fprintf(ficresstdeij," %9.4f (%.4f)", eip, sqrt(vip));
      }
      fprintf(ficresstdeij,"\n");
      fprintf(ficrescveij,"%3.0f",age );
      for(i=1; i<=nlstate;i++)
        for(j=1; j<=nlstate;j++){
-         cptj= (j-1)*nlstate+i;
+                                 cptj= (j-1)*nlstate+i;
-         for(i2=1; i2<=nlstate;i2++)
+                                 for(i2=1; i2<=nlstate;i2++)
-           for(j2=1; j2<=nlstate;j2++){
+                                         for(j2=1; j2<=nlstate;j2++){
-             cptj2= (j2-1)*nlstate+i2;
+                                                 cptj2= (j2-1)*nlstate+i2;
-             if(cptj2 <= cptj)
+                                                 if(cptj2 <= cptj)
-               fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]);
+                                                         fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]);
-           }
+                                         }
        }
      fprintf(ficrescveij,"\n");
    }
    free_matrix(gm,0,nhstepm,1,nlstate*nlstate);
    free_matrix(gp,0,nhstepm,1,nlstate*nlstate);
- Line 4018  void cvevsij(double ***eij, double x[],
+ Line 4660  void cvevsij(double ***eij, double x[],
    free_ma3x(p3matp,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
    printf("\n");
    fprintf(ficlog,"\n");
    free_vector(xm,1,npar);
    free_vector(xp,1,npar);
    free_matrix(dnewm,1,nlstate*nlstate,1,npar);
    free_matrix(doldm,1,nlstate*nlstate,1,nlstate*nlstate);
    free_ma3x(varhe,1,nlstate*nlstate,1,nlstate*nlstate,(int) bage, (int)fage);
  }
  /************ Variance ******************/
   void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[])
- {
+  {
-   /* Variance of health expectancies */
+    /* Variance of health expectancies */
-   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
+    /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
-   /* double **newm;*/
+    /* double **newm;*/
-   /* int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav)*/
+    /* int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav)*/
-   int movingaverage();
+    /* int movingaverage(); */
-   double **dnewm,**doldm;
+    double **dnewm,**doldm;
-   double **dnewmp,**doldmp;
+    double **dnewmp,**doldmp;
-   int i, j, nhstepm, hstepm, h, nstepm ;
+    int i, j, nhstepm, hstepm, h, nstepm ;
-   int k;
+    int k;
-   double *xp;
+    double *xp;
-   double **gp, **gm;  /* for var eij */
+    double **gp, **gm;  /* for var eij */
-   double ***gradg, ***trgradg; /*for var eij */
+    double ***gradg, ***trgradg; /*for var eij */
-   double **gradgp, **trgradgp; /* for var p point j */
+    double **gradgp, **trgradgp; /* for var p point j */
-   double *gpp, *gmp; /* for var p point j */
+    double *gpp, *gmp; /* for var p point j */
-   double **varppt; /* for var p point j nlstate to nlstate+ndeath */
+    double **varppt; /* for var p point j nlstate to nlstate+ndeath */
-   double ***p3mat;
+    double ***p3mat;
-   double age,agelim, hf;
+    double age,agelim, hf;
-   double ***mobaverage;
+    /* double ***mobaverage; */
-   int theta;
+    int theta;
-   char digit[4];
+    char digit[4];
-   char digitp[25];
+    char digitp[25];
-   char fileresprobmorprev[FILENAMELENGTH];
+    char fileresprobmorprev[FILENAMELENGTH];
-   if(popbased==1){
+    if(popbased==1){
-     if(mobilav!=0)
+      if(mobilav!=0)
-       strcpy(digitp,"-POPULBASED-MOBILAV_");
+        strcpy(digitp,"-POPULBASED-MOBILAV_");
-     else strcpy(digitp,"-POPULBASED-NOMOBIL_");
+      else strcpy(digitp,"-POPULBASED-NOMOBIL_");
-   }
+    }
-   else
+    else
-     strcpy(digitp,"-STABLBASED_");
+      strcpy(digitp,"-STABLBASED_");
-   if (mobilav!=0) {
-     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
-     if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
-       fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
-       printf(" Error in movingaverage mobilav=%d\n",mobilav);
-     }
-   }
-   strcpy(fileresprobmorprev,"PRMORPREV-");
-   sprintf(digit,"%-d",ij);
-   /*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/
-   strcat(fileresprobmorprev,digit); /* Tvar to be done */
-   strcat(fileresprobmorprev,digitp); /* Popbased or not, mobilav or not */
-   strcat(fileresprobmorprev,fileresu);
-   if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) {
-     printf("Problem with resultfile: %s\n", fileresprobmorprev);
-     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
-   }
-   printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
-   fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
-   pstamp(ficresprobmorprev);
-   fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
-   fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
-   for(j=nlstate+1; j<=(nlstate+ndeath);j++){
-     fprintf(ficresprobmorprev," p.%-d SE",j);
-     for(i=1; i<=nlstate;i++)
-       fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j);
-   }
-   fprintf(ficresprobmorprev,"\n");
-   fprintf(ficgp,"\n# Routine varevsij");
-   fprintf(ficgp,"\nunset title \n");
- /* fprintf(fichtm, "#Local time at start: %s", strstart);*/
-   fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
-   fprintf(fichtm,"\n<br>%s  <br>\n",digitp);
- /*   } */
-   varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
-   pstamp(ficresvij);
-   fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are ");
-   if(popbased==1)
-     fprintf(ficresvij,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d\n",mobilav);
-   else
-     fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n");
-   fprintf(ficresvij,"# Age");
-   for(i=1; i<=nlstate;i++)
-     for(j=1; j<=nlstate;j++)
-       fprintf(ficresvij," Cov(e.%1d, e.%1d)",i,j);
-   fprintf(ficresvij,"\n");
-   xp=vector(1,npar);
-   dnewm=matrix(1,nlstate,1,npar);
-   doldm=matrix(1,nlstate,1,nlstate);
-   dnewmp= matrix(nlstate+1,nlstate+ndeath,1,npar);
-   doldmp= matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
-   gradgp=matrix(1,npar,nlstate+1,nlstate+ndeath);
-   gpp=vector(nlstate+1,nlstate+ndeath);
-   gmp=vector(nlstate+1,nlstate+ndeath);
-   trgradgp =matrix(nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
-   if(estepm < stepm){
-     printf ("Problem %d lower than %d\n",estepm, stepm);
-   }
-   else  hstepm=estepm;
-   /* For example we decided to compute the life expectancy with the smallest unit */
-   /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.
-      nhstepm is the number of hstepm from age to agelim
-      nstepm is the number of stepm from age to agelim.
-      Look at function hpijx to understand why because of memory size limitations,
-      we decided (b) to get a life expectancy respecting the most precise curvature of the
-      survival function given by stepm (the optimization length). Unfortunately it
-      means that if the survival funtion is printed every two years of age and if
-      you sum them up and add 1 year (area under the trapezoids) you won't get the same
-      results. So we changed our mind and took the option of the best precision.
-   */
-   hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */
-   agelim = AGESUP;
-   for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
-     nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
-     nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
-     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
-     gradg=ma3x(0,nhstepm,1,npar,1,nlstate);
-     gp=matrix(0,nhstepm,1,nlstate);
-     gm=matrix(0,nhstepm,1,nlstate);
-     for(theta=1; theta <=npar; theta++){
-       for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/
-         xp[i] = x[i] + (i==theta ?delti[theta]:0);
-       }
-       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
-       if (popbased==1) {
-         if(mobilav ==0){
-           for(i=1; i<=nlstate;i++)
-             prlim[i][i]=probs[(int)age][i][ij];
-         }else{ /* mobilav */
-           for(i=1; i<=nlstate;i++)
-             prlim[i][i]=mobaverage[(int)age][i][ij];
-         }
-       }
-       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  /* Returns p3mat[i][j][h] for h=1 to nhstepm */
-       for(j=1; j<= nlstate; j++){
-         for(h=0; h<=nhstepm; h++){
-           for(i=1, gp[h][j]=0.;i<=nlstate;i++)
-             gp[h][j] += prlim[i][i]*p3mat[i][j][h];
-         }
-       }
-       /* Next for computing probability of death (h=1 means
-          computed over hstepm matrices product = hstepm*stepm months)
-          as a weighted average of prlim.
-       */
-       for(j=nlstate+1;j<=nlstate+ndeath;j++){
-         for(i=1,gpp[j]=0.; i<= nlstate; i++)
-           gpp[j] += prlim[i][i]*p3mat[i][j][1];
-       }
-       /* end probability of death */
-       for(i=1; i<=npar; i++) /* Computes gradient x - delta */
-         xp[i] = x[i] - (i==theta ?delti[theta]:0);
-       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij);
-       if (popbased==1) {
-         if(mobilav ==0){
-           for(i=1; i<=nlstate;i++)
-             prlim[i][i]=probs[(int)age][i][ij];
-         }else{ /* mobilav */
-           for(i=1; i<=nlstate;i++)
-             prlim[i][i]=mobaverage[(int)age][i][ij];
-         }
-       }
-       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
-       for(j=1; j<= nlstate; j++){  /* Sum of wi * eij = e.j */
-         for(h=0; h<=nhstepm; h++){
-           for(i=1, gm[h][j]=0.;i<=nlstate;i++)
-             gm[h][j] += prlim[i][i]*p3mat[i][j][h];
-         }
-       }
-       /* This for computing probability of death (h=1 means
-          computed over hstepm matrices product = hstepm*stepm months)
-          as a weighted average of prlim.
-       */
-       for(j=nlstate+1;j<=nlstate+ndeath;j++){
-         for(i=1,gmp[j]=0.; i<= nlstate; i++)
-          gmp[j] += prlim[i][i]*p3mat[i][j][1];
-       }
-       /* end probability of death */
-       for(j=1; j<= nlstate; j++) /* vareij */
-         for(h=0; h<=nhstepm; h++){
-           gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
-         }
-       for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */
-         gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
-       }
-     } /* End theta */
-     trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */
-     for(h=0; h<=nhstepm; h++) /* veij */
-       for(j=1; j<=nlstate;j++)
-         for(theta=1; theta <=npar; theta++)
-           trgradg[h][j][theta]=gradg[h][theta][j];
-     for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */
-       for(theta=1; theta <=npar; theta++)
-         trgradgp[j][theta]=gradgp[theta][j];
-     hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
-     for(i=1;i<=nlstate;i++)
-       for(j=1;j<=nlstate;j++)
-         vareij[i][j][(int)age] =0.;
-     for(h=0;h<=nhstepm;h++){
+    /* if (mobilav!=0) { */
-       for(k=0;k<=nhstepm;k++){
+    /*   mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
-         matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);
+    /*   if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ */
-         matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);
+    /*     fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); */
-         for(i=1;i<=nlstate;i++)
+    /*     printf(" Error in movingaverage mobilav=%d\n",mobilav); */
-           for(j=1;j<=nlstate;j++)
+    /*   } */
-             vareij[i][j][(int)age] += doldm[i][j]*hf*hf;
+    /* } */
-       }
-     }
+    strcpy(fileresprobmorprev,"PRMORPREV-");
+    sprintf(digit,"%-d",ij);
+    /*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/
+    strcat(fileresprobmorprev,digit); /* Tvar to be done */
+    strcat(fileresprobmorprev,digitp); /* Popbased or not, mobilav or not */
+    strcat(fileresprobmorprev,fileresu);
+    if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) {
+      printf("Problem with resultfile: %s\n", fileresprobmorprev);
+      fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
+    }
+    printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
+    fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
+    pstamp(ficresprobmorprev);
+    fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
+    fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
+    for(j=nlstate+1; j<=(nlstate+ndeath);j++){
+      fprintf(ficresprobmorprev," p.%-d SE",j);
+      for(i=1; i<=nlstate;i++)
+        fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j);
+    }
+    fprintf(ficresprobmorprev,"\n");
+    fprintf(ficgp,"\n# Routine varevsij");
+    fprintf(ficgp,"\nunset title \n");
+    /* fprintf(fichtm, "#Local time at start: %s", strstart);*/
+    fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
+    fprintf(fichtm,"\n<br>%s  <br>\n",digitp);
+    /*   } */
+    varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
+    pstamp(ficresvij);
+    fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are ");
+    if(popbased==1)
+      fprintf(ficresvij,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d\n",mobilav);
+    else
+      fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n");
+    fprintf(ficresvij,"# Age");
+    for(i=1; i<=nlstate;i++)
+      for(j=1; j<=nlstate;j++)
+        fprintf(ficresvij," Cov(e.%1d, e.%1d)",i,j);
+    fprintf(ficresvij,"\n");
+    xp=vector(1,npar);
+    dnewm=matrix(1,nlstate,1,npar);
+    doldm=matrix(1,nlstate,1,nlstate);
+    dnewmp= matrix(nlstate+1,nlstate+ndeath,1,npar);
+    doldmp= matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
+    gradgp=matrix(1,npar,nlstate+1,nlstate+ndeath);
+    gpp=vector(nlstate+1,nlstate+ndeath);
+    gmp=vector(nlstate+1,nlstate+ndeath);
+    trgradgp =matrix(nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
-     /* pptj */
+    if(estepm < stepm){
-     matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);
+      printf ("Problem %d lower than %d\n",estepm, stepm);
-     matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);
+    }
-     for(j=nlstate+1;j<=nlstate+ndeath;j++)
+    else  hstepm=estepm;
-       for(i=nlstate+1;i<=nlstate+ndeath;i++)
+    /* For example we decided to compute the life expectancy with the smallest unit */
-         varppt[j][i]=doldmp[j][i];
+    /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.
-     /* end ppptj */
+       nhstepm is the number of hstepm from age to agelim
-     /*  x centered again */
+       nstepm is the number of stepm from age to agelim.
+       Look at function hpijx to understand why because of memory size limitations,
-     prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij);
+       we decided (b) to get a life expectancy respecting the most precise curvature of the
+       survival function given by stepm (the optimization length). Unfortunately it
-     if (popbased==1) {
+       means that if the survival funtion is printed every two years of age and if
-       if(mobilav ==0){
+       you sum them up and add 1 year (area under the trapezoids) you won't get the same
-         for(i=1; i<=nlstate;i++)
+       results. So we changed our mind and took the option of the best precision.
-           prlim[i][i]=probs[(int)age][i][ij];
+    */
-       }else{ /* mobilav */
+    hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */
-         for(i=1; i<=nlstate;i++)
+    agelim = AGESUP;
-           prlim[i][i]=mobaverage[(int)age][i][ij];
+    for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
-       }
+      nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
-     }
+      nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
+      p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
-     /* This for computing probability of death (h=1 means
+      gradg=ma3x(0,nhstepm,1,npar,1,nlstate);
-        computed over hstepm (estepm) matrices product = hstepm*stepm months)
+      gp=matrix(0,nhstepm,1,nlstate);
-        as a weighted average of prlim.
+      gm=matrix(0,nhstepm,1,nlstate);
+      for(theta=1; theta <=npar; theta++){
+        for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/
+          xp[i] = x[i] + (i==theta ?delti[theta]:0);
+        }
+        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
+        if (popbased==1) {
+          if(mobilav ==0){
+            for(i=1; i<=nlstate;i++)
+              prlim[i][i]=probs[(int)age][i][ij];
+          }else{ /* mobilav */
+            for(i=1; i<=nlstate;i++)
+              prlim[i][i]=mobaverage[(int)age][i][ij];
+          }
+        }
+        hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  /* Returns p3mat[i][j][h] for h=1 to nhstepm */
+        for(j=1; j<= nlstate; j++){
+          for(h=0; h<=nhstepm; h++){
+            for(i=1, gp[h][j]=0.;i<=nlstate;i++)
+              gp[h][j] += prlim[i][i]*p3mat[i][j][h];
+          }
+        }
+        /* Next for computing probability of death (h=1 means
+           computed over hstepm matrices product = hstepm*stepm months)
+           as a weighted average of prlim.
+        */
+        for(j=nlstate+1;j<=nlstate+ndeath;j++){
+          for(i=1,gpp[j]=0.; i<= nlstate; i++)
+            gpp[j] += prlim[i][i]*p3mat[i][j][1];
+        }
+        /* end probability of death */
+        for(i=1; i<=npar; i++) /* Computes gradient x - delta */
+          xp[i] = x[i] - (i==theta ?delti[theta]:0);
+        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij);
+        if (popbased==1) {
+          if(mobilav ==0){
+            for(i=1; i<=nlstate;i++)
+              prlim[i][i]=probs[(int)age][i][ij];
+          }else{ /* mobilav */
+            for(i=1; i<=nlstate;i++)
+              prlim[i][i]=mobaverage[(int)age][i][ij];
+          }
+        }
+        hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
+        for(j=1; j<= nlstate; j++){  /* Sum of wi * eij = e.j */
+          for(h=0; h<=nhstepm; h++){
+            for(i=1, gm[h][j]=0.;i<=nlstate;i++)
+              gm[h][j] += prlim[i][i]*p3mat[i][j][h];
+          }
+        }
+        /* This for computing probability of death (h=1 means
+           computed over hstepm matrices product = hstepm*stepm months)
+           as a weighted average of prlim.
+        */
+        for(j=nlstate+1;j<=nlstate+ndeath;j++){
+          for(i=1,gmp[j]=0.; i<= nlstate; i++)
+            gmp[j] += prlim[i][i]*p3mat[i][j][1];
+        }
+        /* end probability of death */
+        for(j=1; j<= nlstate; j++) /* vareij */
+          for(h=0; h<=nhstepm; h++){
+            gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
+          }
+        for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */
+          gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
+        }
+      } /* End theta */
+      trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */
+      for(h=0; h<=nhstepm; h++) /* veij */
+        for(j=1; j<=nlstate;j++)
+          for(theta=1; theta <=npar; theta++)
+            trgradg[h][j][theta]=gradg[h][theta][j];
+      for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */
+        for(theta=1; theta <=npar; theta++)
+          trgradgp[j][theta]=gradgp[theta][j];
+      hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
+      for(i=1;i<=nlstate;i++)
+        for(j=1;j<=nlstate;j++)
+          vareij[i][j][(int)age] =0.;
+      for(h=0;h<=nhstepm;h++){
+        for(k=0;k<=nhstepm;k++){
+          matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);
+          matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);
+          for(i=1;i<=nlstate;i++)
+            for(j=1;j<=nlstate;j++)
+              vareij[i][j][(int)age] += doldm[i][j]*hf*hf;
+        }
+      }
+      /* pptj */
+      matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);
+      matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);
+      for(j=nlstate+1;j<=nlstate+ndeath;j++)
+        for(i=nlstate+1;i<=nlstate+ndeath;i++)
+          varppt[j][i]=doldmp[j][i];
+      /* end ppptj */
+      /*  x centered again */
+      prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij);
+      if (popbased==1) {
+        if(mobilav ==0){
+          for(i=1; i<=nlstate;i++)
+            prlim[i][i]=probs[(int)age][i][ij];
+        }else{ /* mobilav */
+          for(i=1; i<=nlstate;i++)
+            prlim[i][i]=mobaverage[(int)age][i][ij];
+        }
+      }
+      /* This for computing probability of death (h=1 means
+         computed over hstepm (estepm) matrices product = hstepm*stepm months)
+         as a weighted average of prlim.
+      */
+      hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);
+      for(j=nlstate+1;j<=nlstate+ndeath;j++){
+        for(i=1,gmp[j]=0.;i<= nlstate; i++)
+          gmp[j] += prlim[i][i]*p3mat[i][j][1];
+      }
+      /* end probability of death */
+      fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);
+      for(j=nlstate+1; j<=(nlstate+ndeath);j++){
+        fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));
+        for(i=1; i<=nlstate;i++){
+          fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]);
+        }
+      }
+      fprintf(ficresprobmorprev,"\n");
+      fprintf(ficresvij,"%.0f ",age );
+      for(i=1; i<=nlstate;i++)
+        for(j=1; j<=nlstate;j++){
+          fprintf(ficresvij," %.4f", vareij[i][j][(int)age]);
+        }
+      fprintf(ficresvij,"\n");
+      free_matrix(gp,0,nhstepm,1,nlstate);
+      free_matrix(gm,0,nhstepm,1,nlstate);
+      free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate);
+      free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);
+      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+    } /* End age */
+    free_vector(gpp,nlstate+1,nlstate+ndeath);
+    free_vector(gmp,nlstate+1,nlstate+ndeath);
+    free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath);
+    free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
+    /* fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); */
+    fprintf(ficgp,"\nunset parametric;unset label; set ter svg size 640, 480");
+    /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
+    fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
+    fprintf(ficgp,"\nset out \"%s%s.svg\";",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
+    /*   fprintf(ficgp,"\n plot \"%s\"  u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */
+    /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
+    /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
+    fprintf(ficgp,"\n plot \"%s\"  u 1:($3) not w l lt 1 ",subdirf(fileresprobmorprev));
+    fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev));
+    fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev));
+    fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));
+    fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.svg\"> <br>\n", estepm,subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
+    /*  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.svg\"> <br>\n", stepm,YEARM,digitp,digit);
      */
-     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);
+    /*   fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.svg\";replot;",digitp,optionfilefiname,digit); */
-     for(j=nlstate+1;j<=nlstate+ndeath;j++){
+    fprintf(ficgp,"\nset out;\nset out \"%s%s.svg\";replot;set out;\n",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
-       for(i=1,gmp[j]=0.;i<= nlstate; i++)
-         gmp[j] += prlim[i][i]*p3mat[i][j][1];
-     }
-     /* end probability of death */
-     fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);
-     for(j=nlstate+1; j<=(nlstate+ndeath);j++){
-       fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));
-       for(i=1; i<=nlstate;i++){
-         fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]);
-       }
-     }
-     fprintf(ficresprobmorprev,"\n");
-     fprintf(ficresvij,"%.0f ",age );
+    free_vector(xp,1,npar);
-     for(i=1; i<=nlstate;i++)
+    free_matrix(doldm,1,nlstate,1,nlstate);
-       for(j=1; j<=nlstate;j++){
+    free_matrix(dnewm,1,nlstate,1,npar);
-         fprintf(ficresvij," %.4f", vareij[i][j][(int)age]);
+    free_matrix(doldmp,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
-       }
+    free_matrix(dnewmp,nlstate+1,nlstate+ndeath,1,npar);
-     fprintf(ficresvij,"\n");
+    free_matrix(varppt,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
-     free_matrix(gp,0,nhstepm,1,nlstate);
+    /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
-     free_matrix(gm,0,nhstepm,1,nlstate);
+    fclose(ficresprobmorprev);
-     free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate);
+    fflush(ficgp);
-     free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);
+    fflush(fichtm);
-     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+  }  /* end varevsij */
-   } /* End age */
-   free_vector(gpp,nlstate+1,nlstate+ndeath);
-   free_vector(gmp,nlstate+1,nlstate+ndeath);
-   free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath);
-   free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
-   /* fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); */
-   fprintf(ficgp,"\nunset parametric;unset label; set ter svg size 640, 480");
-   /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
-   fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
-   fprintf(ficgp,"\nset out \"%s%s.svg\";",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
- /*   fprintf(ficgp,"\n plot \"%s\"  u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */
- /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
- /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
-   fprintf(ficgp,"\n plot \"%s\"  u 1:($3) not w l lt 1 ",subdirf(fileresprobmorprev));
-   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev));
-   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev));
-   fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));
-   fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.svg\"> <br>\n", estepm,subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
-   /*  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.svg\"> <br>\n", stepm,YEARM,digitp,digit);
- */
- /*   fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.svg\";replot;",digitp,optionfilefiname,digit); */
-   fprintf(ficgp,"\nset out;\nset out \"%s%s.svg\";replot;set out;\n",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
-   free_vector(xp,1,npar);
-   free_matrix(doldm,1,nlstate,1,nlstate);
-   free_matrix(dnewm,1,nlstate,1,npar);
-   free_matrix(doldmp,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
-   free_matrix(dnewmp,nlstate+1,nlstate+ndeath,1,npar);
-   free_matrix(varppt,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
-   if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
-   fclose(ficresprobmorprev);
-   fflush(ficgp);
-   fflush(fichtm);
- }  /* end varevsij */
  /************ Variance of prevlim ******************/
   void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[])
- Line 4800  To be simple, these graphs help to under
+ Line 5442  To be simple, these graphs help to under
  void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \
                    int lastpass, int stepm, int weightopt, char model[],\
                    int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
-                   int popforecast, int estepm ,\
+                   int popforecast, int prevfcast, int backcast, int estepm , \
-                   double jprev1, double mprev1,double anprev1, \
+                   double jprev1, double mprev1,double anprev1, double dateprev1, \
-                   double jprev2, double mprev2,double anprev2){
+                   double jprev2, double mprev2,double anprev2, double dateprev2){
    int jj1, k1, i1, cpt;
     fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
     <li><a href='#secondorder'>Result files (second order (variance)</a>\n \
  </ul>");
-    fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n \
+    fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n");
-  - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ",
+    fprintf(fichtm,"<li>- Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file)<br/>\n",
-            jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_"));
+            jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm"));
+    fprintf(fichtm,"<li> - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file) ",
+            jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTM_",".htm"),subdirfext3(optionfilefiname,"PHTM_",".htm"));
+    fprintf(fichtm,",  <a href=\"%s\">%s</a> (text file) <br>\n",subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_"));
     fprintf(fichtm,"\
   - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
             stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_"));
     fprintf(fichtm,"\
+  - Estimated back transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
+            stepm,subdirf2(fileresu,"PIJB_"),subdirf2(fileresu,"PIJB_"));
+    fprintf(fichtm,"\
   - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
             subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_"));
     fprintf(fichtm,"\
-  - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
+  - Period (stable) back prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
+            subdirf2(fileresu,"PLB_"),subdirf2(fileresu,"PLB_"));
+    fprintf(fichtm,"\
+  - (a) Life expectancies by health status at initial age, e<sub>i.</sub> (b) health expectancies by health status at initial age, e<sub>ij</sub> . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
     <a href=\"%s\">%s</a> <br>\n",
             estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_"));
-    fprintf(fichtm,"\
+    if(prevfcast==1){
-  - Population projections by age and states: \
+      fprintf(fichtm,"\
+  - Prevalence projections by age and states:                            \
     <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_"));
+    }
  fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
- Line 4843  fprintf(fichtm," \n<ul><li><b>Graphs</b>
+ Line 5496  fprintf(fichtm," \n<ul><li><b>Graphs</b>
         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
       }
       /* aij, bij */
-      fprintf(fichtm,"<br>- Logit model, for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
+      fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
- <img src=\"%s_%d-1.svg\">",subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
+ <img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
       /* Pij */
-      fprintf(fichtm,"<br>\n- Pij or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
+      fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
  <img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
       /* Quasi-incidences */
-      fprintf(fichtm,"<br>\n- Iij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
+      fprintf(fichtm,"<br>\n- I<sub>ij</sub> or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
   before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too,\
-  incidence (rates) are the limit when h tends to zero of the ratio of the probability hPij \
+  incidence (rates) are the limit when h tends to zero of the ratio of the probability  <sub>h</sub>P<sub>ij</sub> \
- divided by h: hPij/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \
+ divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \
  <img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
       /* Survival functions (period) in state j */
       for(cpt=1; cpt<=nlstate;cpt++){
- Line 4867  divided by h: hPij/h : <a href=\"%s_%d-3
+ Line 5520  divided by h: hPij/h : <a href=\"%s_%d-3
       }
       /* Period (stable) prevalence in each health state */
       for(cpt=1; cpt<=nlstate;cpt++){
-        fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
+        fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
  <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);
       }
+     if(backcast==1){
+      /* Period (stable) back prevalence in each health state */
+      for(cpt=1; cpt<=nlstate;cpt++){
+        fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
+ <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1);
+      }
+     }
+     if(prevfcast==1){
+       /* Projection of prevalence up to period (stable) prevalence in each health state */
+       for(cpt=1; cpt<=nlstate;cpt++){
+         fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
+ <img src=\"%s_%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1);
+       }
+     }
       for(cpt=1; cpt<=nlstate;cpt++) {
         fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \
  <img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1);
- Line 4940  See page 'Matrix of variance-covariance
+ Line 5608  See page 'Matrix of variance-covariance
         fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
       }
       for(cpt=1; cpt<=nlstate;cpt++) {
-        fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \
+        fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \
- prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d%d.svg\"> %s_%d-%d.svg <br>\
+ prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d.svg\"> %s_%d-%d.svg</a>\n <br>\
  <img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1);
       }
       fprintf(fichtm,"\n<br>- Total life expectancy by age and \
  health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
  true period expectancies (those weighted with period prevalences are also\
   drawn in addition to the population based expectancies computed using\
-  observed and cahotic prevalences:  <a href=\"%s_%d.svg\">%s_%d.svg<br>\
+  observed and cahotic prevalences:  <a href=\"%s_%d.svg\">%s_%d.svg</a>\n<br>\
  <img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1);
     /* } /\* end i1 *\/ */
   }/* End k1 */
- Line 4957  true period expectancies (those weighted
+ Line 5625  true period expectancies (those weighted
  }
  /******************* Gnuplot file **************/
- void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
+  void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){
    char dirfileres[132],optfileres[132];
    int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
+   int lv=0, vlv=0, kl=0;
    int ng=0;
    int vpopbased;
  /*   if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
- Line 4984  void printinggnuplot(char fileresu[], ch
+ Line 5653  void printinggnuplot(char fileresu[], ch
  /* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)"  */
      /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */
      fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_"));
-     fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$12):5 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk));
+     fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$13):6 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk));
      fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_"));
-     fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$12):4 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk));
+     fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$13):5 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk));
      for (i=1; i<= nlstate ; i ++) {
        fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i);
        fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot  \"%s\"",subdirf(fileresilk));
-       fprintf(ficgp,"  u  2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);
+       fprintf(ficgp,"  u  2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);
        for (j=2; j<= nlstate+ndeath ; j ++) {
-         fprintf(ficgp,",\\\n \"\" u  2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
+         fprintf(ficgp,",\\\n \"\" u  2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
        }
        fprintf(ficgp,";\nset out; unset ylabel;\n");
      }
- Line 5005  void printinggnuplot(char fileresu[], ch
+ Line 5674  void printinggnuplot(char fileresu[], ch
    strcpy(dirfileres,optionfilefiname);
    strcpy(optfileres,"vpl");
   /* 1eme*/
-   fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files\n");
+   for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */
-   for (cpt=1; cpt<= nlstate ; cpt ++) {
+     for (k1=1; k1<= m ; k1 ++) { /* For each combination of covariate */
-     for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
+       /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
+       fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files ");
+       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+         lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+         /* 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[lv]][lv];
+         fprintf(ficgp," V%d=%d ",k,vlv);
+       }
+       fprintf(ficgp,"\n#\n");
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
       fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \n\
- Line 5029  plot [%.f:%.f] \"%s\" every :::%d::%d u
+ Line 5709  plot [%.f:%.f] \"%s\" every :::%d::%d u
         if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
         else fprintf(ficgp," %%*lf (%%*lf)");
       }
-      fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
+      fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
+                  if(backcast==1){
+                          fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt);
+                  }
       fprintf(ficgp,"\nset out \n");
      } /* k1 */
    } /* cpt */
    /*2 eme*/
-   fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n");
    for (k1=1; k1<= m ; k1 ++) {
+       fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
+       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+         lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+         /* 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[lv]][lv];
+         fprintf(ficgp," V%d=%d ",k,vlv);
+       }
+       fprintf(ficgp,"\n#\n");
      fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
      for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
        if(vpopbased==0)
- Line 5068  plot [%.f:%.f] \"%s\" every :::%d::%d u
+ Line 5761  plot [%.f:%.f] \"%s\" every :::%d::%d u
      } /* vpopbased */
      fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
    } /* k1 */
    /*3eme*/
    for (k1=1; k1<= m ; k1 ++) {
      for (cpt=1; cpt<= nlstate ; cpt ++) {
+       fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files:  cov=%d state=%d",k1, cpt);
+       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+         lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+         /* 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[lv]][lv];
+         fprintf(ficgp," V%d=%d ",k,vlv);
+       }
+       fprintf(ficgp,"\n#\n");
        /*       k=2+nlstate*(2*cpt-2); */
        k=2+(nlstate+1)*(cpt-1);
        fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
- Line 5097  plot [%.f:%.f] \"%s\" every :::%d::%d u
+ Line 5802  plot [%.f:%.f] \"%s\" every :::%d::%d u
    /* Survival functions (period) from state i in state j by initial state i */
    for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
-       k=3;
+       fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
-       fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'lij' files, cov=%d state=%d",k1, cpt);
+       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+         lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+         /* 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[lv]][lv];
+         fprintf(ficgp," V%d=%d ",k,vlv);
+       }
+       fprintf(ficgp,"\n#\n");
        fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
  set ter svg size 640, 480\n\
  unset log y\n\
  plot [%.f:%.f]  ", ageminpar, agemaxpar);
+       k=3;
        for (i=1; i<= nlstate ; i ++){
-         if(i==1)
+         if(i==1){
            fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
-         else
+         }else{
            fprintf(ficgp,", '' ");
+         }
          l=(nlstate+ndeath)*(i-1)+1;
          fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
          for (j=2; j<= nlstate+ndeath ; j ++)
- Line 5122  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 5838  plot [%.f:%.f]  ", ageminpar, agemaxpar)
    /* Survival functions (period) from state i in state j by final state j */
    for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
-       k=3;
        fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
+       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+         lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+         /* 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[lv]][lv];
+         fprintf(ficgp," V%d=%d ",k,vlv);
+       }
+       fprintf(ficgp,"\n#\n");
        fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
  set ter svg size 640, 480\n\
  unset log y\n\
  plot [%.f:%.f]  ", ageminpar, agemaxpar);
+       k=3;
        for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
          if(j==1)
            fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
- Line 5154  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 5880  plot [%.f:%.f]  ", ageminpar, agemaxpar)
    } /* end covariate */
    /* CV preval stable (period) for each covariate */
-   for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
+   for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
-       k=3;
+       fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
-       fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt);
+       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+         lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+         /* 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[lv]][lv];
+         fprintf(ficgp," V%d=%d ",k,vlv);
+       }
+       fprintf(ficgp,"\n#\n");
        fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);
        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
  set ter svg size 640, 480\n\
  unset log y\n\
  plot [%.f:%.f]  ", ageminpar, agemaxpar);
+       k=3; /* Offset */
        for (i=1; i<= nlstate ; i ++){
          if(i==1)
            fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
- Line 5177  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 5913  plot [%.f:%.f]  ", ageminpar, agemaxpar)
        fprintf(ficgp,"\nset out\n");
      } /* end cpt state*/
    } /* end covariate */
+   if(backcast == 1){
+     /* CV back preval stable (period) for each covariate */
+     for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
+       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
+         fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
+         for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
+           lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+           /* 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[lv]][lv];
+           fprintf(ficgp," V%d=%d ",k,vlv);
+         }
+         fprintf(ficgp,"\n#\n");
+         fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
+         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
+ set ter svg size 640, 480\n                                             \
+ unset log y\n                                                           \
+ plot [%.f:%.f]  ", ageminpar, agemaxpar);
+         k=3; /* Offset */
+         for (i=1; i<= nlstate ; i ++){
+           if(i==1)
+             fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));
+           else
+             fprintf(ficgp,", '' ");
+           /* l=(nlstate+ndeath)*(i-1)+1; */
+           l=(nlstate+ndeath)*(cpt-1)+1;
+           /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */
+           /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */
+           fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+(cpt-1)+i-1); /* a vérifier */
+           /* for (j=2; j<= nlstate ; j ++) */
+           /*    fprintf(ficgp,"+$%d",k+l+j-1); */
+           /*    /\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */
+           fprintf(ficgp,") t \"bprev(%d,%d)\" w l",i,cpt);
+         } /* nlstate */
+         fprintf(ficgp,"\nset out\n");
+       } /* end cpt state*/
+     } /* end covariate */
+   } /* End if backcast */
+   if(prevfcast==1){
+     /* Projection from cross-sectional to stable (period) for each covariate */
+     for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
+       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
+         fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
+         for (k=1; k<=cptcoveff; k++){    /* For each correspondig covariate value  */
+           lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+           /* 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[lv]][lv];
+           fprintf(ficgp," V%d=%d ",k,vlv);
+         }
+         fprintf(ficgp,"\n#\n");
+         fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
+         fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
+         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
+ set ter svg size 640, 480\n\
+ unset log y\n\
+ plot [%.f:%.f]  ", ageminpar, agemaxpar);
+         for (i=1; i<= nlstate+1 ; i ++){  /* nlstate +1 p11 p21 p.1 */
+           /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+           /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */
+           /*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+           /*#   1       2   3    4    5      6  7   8   9   10   11 12  13   14  15 */
+           if(i==1){
+             fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_"));
+           }else{
+             fprintf(ficgp,",\\\n '' ");
+           }
+           if(cptcoveff ==0){ /* No covariate */
+             fprintf(ficgp," u 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 */
+             if(i==nlstate+1)
+               fprintf(ficgp," $%d/(1.-$%d)) t 'p.%d' with line ", \
++(cpt-1)*(nlstate+1)+1+(i-1),  2+1+(i-1)+(nlstate+1)*nlstate,cpt );
+             else
+               fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \
++(cpt-1)*(nlstate+1)+1+(i-1),  2+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+           }else{
+             fprintf(ficgp,"u 6:(("); /* Age is in 6 */
+             /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+             /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */
+             kl=0;
+             for (k=1; k<=cptcoveff; k++){    /* For each covariate  */
+               lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+               /* 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[lv]][lv];
+               kl++;
+               /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
+               /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */
+               /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */
+               /* ''  u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
+               if(k==cptcoveff)
+                 if(i==nlstate+1)
+                   fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \
++(cpt-1)*(nlstate+1)+1+(i-1),  6+1+(i-1)+(nlstate+1)*nlstate,cpt );
+                 else
+                   fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \
++(cpt-1)*(nlstate+1)+1+(i-1),  6+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+               else{
+                 fprintf(ficgp,"$%d==%d && $%d==%d && ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv]);
+                 kl++;
+               }
+             } /* end covariate */
+           } /* end if covariate */
+         } /* nlstate */
+         fprintf(ficgp,"\nset out\n");
+       } /* end cpt state*/
+     } /* end covariate */
+   } /* End if prevfcast */
    /* proba elementaires */
    fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
- Line 5267  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6121  plot [%.f:%.f]  ", ageminpar, agemaxpar)
                 else
                   fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
               }
-              if(ng != 1){
+            }else{
-                fprintf(ficgp,")/(1");
+              i=i-ncovmodel;
+              if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
+                fprintf(ficgp," (1.");
+            }
+            if(ng != 1){
+              fprintf(ficgp,")/(1");
-                for(k1=1; k1 <=nlstate; k1++){
+              for(k1=1; k1 <=nlstate; k1++){
-                  if(nagesqr==0)
+                if(nagesqr==0)
-                    fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+                  fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
-                  else /* nagesqr =1 */
+                else /* nagesqr =1 */
-                    fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
+                  fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
-                  ij=1;
+                ij=1;
-                  for(j=3; j <=ncovmodel-nagesqr; j++){
+                for(j=3; j <=ncovmodel-nagesqr; j++){
-                    if(ij <=cptcovage) { /* Bug valgrind */
+                  if(ij <=cptcovage) { /* Bug valgrind */
-                      if((j-2)==Tage[ij]) { /* Bug valgrind */
+                    if((j-2)==Tage[ij]) { /* Bug valgrind */
-                        fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+                      fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
-                        /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+                      /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
-                        ij++;
+                      ij++;
-                      }
                     }
-                    else
-                      fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
                   }
-                  fprintf(ficgp,")");
+                  else
+                    fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
                 }
                 fprintf(ficgp,")");
-                if(ng ==2)
-                  fprintf(ficgp," t \"p%d%d\" ", k2,k);
-                else /* ng= 3 */
-                  fprintf(ficgp," t \"i%d%d\" ", k2,k);
-              }else{ /* end ng <> 1 */
-                fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
               }
-              if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
+              fprintf(ficgp,")");
-              i=i+ncovmodel;
+              if(ng ==2)
+                fprintf(ficgp," t \"p%d%d\" ", k2,k);
+              else /* ng= 3 */
+                fprintf(ficgp," t \"i%d%d\" ", k2,k);
+            }else{ /* end ng <> 1 */
+              if( k !=k2) /* logit p11 is hard to draw */
+                fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
             }
+            if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)
+              fprintf(ficgp,",");
+            if (ng == 1 && k!=k2 && (k+k2)!= (nlstate*2+ndeath))
+              fprintf(ficgp,",");
+            i=i+ncovmodel;
           } /* end k */
         } /* end k2 */
         fprintf(ficgp,"\n set out\n");
- Line 5312  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6175  plot [%.f:%.f]  ", ageminpar, agemaxpar)
  /*************** Moving average **************/
- int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav){
+                 /* int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav, double bageout, double fageout){ */
+ int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav){
    int i, cpt, cptcod;
    int modcovmax =1;
    int mobilavrange, mob;
    double age;
+   int iage=0;
+   double *sumnewp, *sumnewm;
+   double *agemingood, *agemaxgood; /* Currently identical for all covariates */
+   sumnewp = vector(1,modcovmax);
+   sumnewm = vector(1,modcovmax);
+   agemingood = vector(1,modcovmax);
+   agemaxgood = vector(1,modcovmax);
    modcovmax=2*cptcoveff;/* Max number of modalities. We suppose
-                            a covariate has 2 modalities */
+                            a covariate has 2 modalities, should be equal to ncovcombmax  */
    if (cptcovn<1) modcovmax=1; /* At least 1 pass */
    if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){
      if(mobilav==1) mobilavrange=5; /* default */
      else mobilavrange=mobilav;
- Line 5339  int movingaverage(double ***probs, doubl
+ Line 6212  int movingaverage(double ***probs, doubl
          for (i=1; i<=nlstate;i++){
            for (cptcod=1;cptcod<=modcovmax;cptcod++){
              mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod];
-               for (cpt=1;cpt<=(mob-1)/2;cpt++){
+             for (cpt=1;cpt<=(mob-1)/2;cpt++){
-                 mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod];
+               mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod];
-                 mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod];
+               mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod];
-               }
+             }
              mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob;
            }
          }
        }/* end age */
      }/* end mob */
-   }else return -1;
+   }else
+     return -1;
+   for (cptcod=1;cptcod<=modcovmax;cptcod++){
+     /* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */
+     agemingood[cptcod]=fage+(mob-1)/2;
+     for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, finding the youngest wrong */
+       sumnewm[cptcod]=0.;
+       for (i=1; i<=nlstate;i++){
+         sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+       }
+       if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
+         agemingood[cptcod]=age;
+       }else{ /* bad */
+         for (i=1; i<=nlstate;i++){
+           mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
+         } /* i */
+       } /* end bad */
+     }/* age */
+     /* From youngest, finding the oldest wrong */
+     agemaxgood[cptcod]=bage+(mob-1)/2;
+     for (age=bage+(mob-1)/2; age<=fage; age++){
+       sumnewm[cptcod]=0.;
+       for (i=1; i<=nlstate;i++){
+         sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+       }
+       if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
+         agemaxgood[cptcod]=age;
+       }else{ /* bad */
+         for (i=1; i<=nlstate;i++){
+           mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
+         } /* i */
+       } /* end bad */
+     }/* age */
+     for (age=bage; age<=fage; age++){
+       printf("%d %d ", cptcod, (int)age);
+       sumnewp[cptcod]=0.;
+       sumnewm[cptcod]=0.;
+       for (i=1; i<=nlstate;i++){
+         sumnewp[cptcod]+=probs[(int)age][i][cptcod];
+         sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+         printf("%.4f %.4f ",probs[(int)age][i][cptcod], mobaverage[(int)age][i][cptcod]);
+       }
+       printf("%.4f %.4f \n",sumnewp[cptcod], sumnewm[cptcod]);
+     }
+     printf("\n");
+     /* brutal averaging */
+     for (i=1; i<=nlstate;i++){
+       for (age=1; age<=bage; age++){
+         mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
+         printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]);
+       }
+       for (age=fage; age<=AGESUP; age++){
+         mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
+         printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]);
+       }
+     } /* end i status */
+     for (i=nlstate+1; i<=nlstate+ndeath;i++){
+       for (age=1; age<=AGESUP; age++){
+         /*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/
+         mobaverage[(int)age][i][cptcod]=0.;
+       }
+     }
+   }/* end cptcod */
+   free_vector(sumnewm,1, modcovmax);
+   free_vector(sumnewp,1, modcovmax);
+   free_vector(agemaxgood,1, modcovmax);
+   free_vector(agemingood,1, modcovmax);
    return 0;
  }/* End movingaverage */
  /************** Forecasting ******************/
  void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
- Line 5365  void prevforecast(char fileres[], double
+ Line 6304  void prevforecast(char fileres[], double
    double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
    double *popeffectif,*popcount;
    double ***p3mat;
-   double ***mobaverage;
+   /* double ***mobaverage; */
    char fileresf[FILENAMELENGTH];
    agelim=AGESUP;
-   prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+   /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
+      in each health status at the date of interview (if between dateprev1 and dateprev2).
+      We still use firstpass and lastpass as another selection.
+   */
+   /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ */
+   /*          firstpass, lastpass,  stepm,  weightopt, model); */
    strcpy(fileresf,"F_");
    strcat(fileresf,fileresu);
- Line 5377  void prevforecast(char fileres[], double
+ Line 6321  void prevforecast(char fileres[], double
      printf("Problem with forecast resultfile: %s\n", fileresf);
      fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf);
    }
-   printf("Computing forecasting: result on file '%s' \n", fileresf);
+   printf("Computing forecasting: result on file '%s', please wait... \n", fileresf);
-   fprintf(ficlog,"Computing forecasting: result on file '%s' \n", fileresf);
+   fprintf(ficlog,"Computing forecasting: result on file '%s', please wait... \n", fileresf);
    if (cptcoveff==0) ncodemax[cptcoveff]=1;
-   if (mobilav!=0) {
-     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
-     if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
-       fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
-       printf(" Error in movingaverage mobilav=%d\n",mobilav);
-     }
-   }
    stepsize=(int) (stepm+YEARM-1)/YEARM;
    if (stepm<=12) stepsize=1;
- Line 5419  void prevforecast(char fileres[], double
+ Line 6356  void prevforecast(char fileres[], double
    for(cptcov=1, k=0;cptcov<=i1;cptcov++){
      for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
        k=k+1;
-       fprintf(ficresf,"\n#******");
+       fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#");
        for(j=1;j<=cptcoveff;j++) {
-         fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+         fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
        }
-       fprintf(ficresf,"******\n");
+       fprintf(ficresf," yearproj age");
-       fprintf(ficresf,"# Covariate valuofcovar yearproj age");
        for(j=1; j<=nlstate+ndeath;j++){
          for(i=1; i<=nlstate;i++)
            fprintf(ficresf," p%d%d",i,j);
          fprintf(ficresf," p.%d",j);
        }
        for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
          fprintf(ficresf,"\n");
          fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);
          for (agec=fage; agec>=(ageminpar-1); agec--){
            nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
            nhstepm = nhstepm/hstepm;
            p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
            oldm=oldms;savm=savms;
            hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
            for (h=0; h<=nhstepm; h++){
              if (h*hstepm/YEARM*stepm ==yearp) {
- Line 5471  void prevforecast(char fileres[], double
+ Line 6406  void prevforecast(char fileres[], double
      } /* end cptcod */
    } /* end  cptcov */
-   if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
    fclose(ficresf);
+   printf("End of Computing forecasting \n");
+   fprintf(ficlog,"End of Computing forecasting\n");
  }
+ /* /\************** Back Forecasting ******************\/ */
+ /* void prevbackforecast(char fileres[], double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){ */
+ /*   /\* back1, year, month, day of starting backection  */
+ /*      agemin, agemax range of age */
+ /*      dateprev1 dateprev2 range of dates during which prevalence is computed */
+ /*      anback2 year of en of backection (same day and month as back1). */
+ /*   *\/ */
+ /*   int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1; */
+ /*   double agec; /\* generic age *\/ */
+ /*   double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; */
+ /*   double *popeffectif,*popcount; */
+ /*   double ***p3mat; */
+ /*   /\* double ***mobaverage; *\/ */
+ /*   char fileresfb[FILENAMELENGTH]; */
+ /*   agelim=AGESUP; */
+ /*   /\* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people */
+ /*      in each health status at the date of interview (if between dateprev1 and dateprev2). */
+ /*      We still use firstpass and lastpass as another selection. */
+ /*   *\/ */
+ /*   /\* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ *\/ */
+ /*   /\*              firstpass, lastpass,  stepm,  weightopt, model); *\/ */
+ /*   prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
+ /*   strcpy(fileresfb,"FB_");  */
+ /*   strcat(fileresfb,fileresu); */
+ /*   if((ficresfb=fopen(fileresfb,"w"))==NULL) { */
+ /*     printf("Problem with back forecast resultfile: %s\n", fileresfb); */
+ /*     fprintf(ficlog,"Problem with back forecast resultfile: %s\n", fileresfb); */
+ /*   } */
+ /*   printf("Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */
+ /*   fprintf(ficlog,"Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */
+ /*   if (cptcoveff==0) ncodemax[cptcoveff]=1; */
+ /*   /\* if (mobilav!=0) { *\/ */
+ /*   /\*   mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */
+ /*   /\*   if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ *\/ */
+ /*   /\*     fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); *\/ */
+ /*   /\*     printf(" Error in movingaverage mobilav=%d\n",mobilav); *\/ */
+ /*   /\*   } *\/ */
+ /*   /\* } *\/ */
+ /*   stepsize=(int) (stepm+YEARM-1)/YEARM; */
+ /*   if (stepm<=12) stepsize=1; */
+ /*   if(estepm < stepm){ */
+ /*     printf ("Problem %d lower than %d\n",estepm, stepm); */
+ /*   } */
+ /*   else  hstepm=estepm;    */
+ /*   hstepm=hstepm/stepm;  */
+ /*   yp1=modf(dateintmean,&yp);/\* extracts integral of datemean in yp  and */
+ /*                                fractional in yp1 *\/ */
+ /*   anprojmean=yp; */
+ /*   yp2=modf((yp1*12),&yp); */
+ /*   mprojmean=yp; */
+ /*   yp1=modf((yp2*30.5),&yp); */
+ /*   jprojmean=yp; */
+ /*   if(jprojmean==0) jprojmean=1; */
+ /*   if(mprojmean==0) jprojmean=1; */
+ /*   i1=cptcoveff; */
+ /*   if (cptcovn < 1){i1=1;} */
+ /*   fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);  */
+ /*   fprintf(ficresfb,"#****** Routine prevbackforecast **\n"); */
+ /*      /\*           if (h==(int)(YEARM*yearp)){ *\/ */
+ /*   for(cptcov=1, k=0;cptcov<=i1;cptcov++){ */
+ /*     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */
+ /*       k=k+1; */
+ /*       fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.jx is weighted by observed prev \n#"); */
+ /*       for(j=1;j<=cptcoveff;j++) { */
+ /*                              fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */
+ /*       } */
+ /*       fprintf(ficresfb," yearbproj age"); */
+ /*       for(j=1; j<=nlstate+ndeath;j++){  */
+ /*                              for(i=1; i<=nlstate;i++)               */
+ /*           fprintf(ficresfb," p%d%d",i,j); */
+ /*                              fprintf(ficresfb," p.%d",j); */
+ /*       } */
+ /*       for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) {  */
+ /*                              /\* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {  *\/ */
+ /*                              fprintf(ficresfb,"\n"); */
+ /*                              fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);    */
+ /*                              for (agec=fage; agec>=(ageminpar-1); agec--){  */
+ /*                                      nhstepm=(int) rint((agelim-agec)*YEARM/stepm);  */
+ /*                                      nhstepm = nhstepm/hstepm;  */
+ /*                                      p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
+ /*                                      oldm=oldms;savm=savms; */
+ /*                                      hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm,oldm,savm, dnewm, doldm, dsavm, k);       */
+ /*                                      for (h=0; h<=nhstepm; h++){ */
+ /*                                              if (h*hstepm/YEARM*stepm ==yearp) { */
+ /*               fprintf(ficresfb,"\n"); */
+ /*               for(j=1;j<=cptcoveff;j++)  */
+ /*                 fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */
+ /*                                                      fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec+h*hstepm/YEARM*stepm); */
+ /*                                              }  */
+ /*                                              for(j=1; j<=nlstate+ndeath;j++) { */
+ /*                                                      ppij=0.; */
+ /*                                                      for(i=1; i<=nlstate;i++) { */
+ /*                                                              if (mobilav==1)  */
+ /*                                                                      ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod]; */
+ /*                                                              else { */
+ /*                                                                      ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod]; */
+ /*                                                              } */
+ /*                                                              if (h*hstepm/YEARM*stepm== yearp) { */
+ /*                                                                      fprintf(ficresfb," %.3f", p3mat[i][j][h]); */
+ /*                                                              } */
+ /*                                                      } /\* end i *\/ */
+ /*                                                      if (h*hstepm/YEARM*stepm==yearp) { */
+ /*                                                              fprintf(ficresfb," %.3f", ppij); */
+ /*                                                      } */
+ /*                                              }/\* end j *\/ */
+ /*                                      } /\* end h *\/ */
+ /*                                      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
+ /*                              } /\* end agec *\/ */
+ /*       } /\* end yearp *\/ */
+ /*     } /\* end cptcod *\/ */
+ /*   } /\* end  cptcov *\/ */
+ /*   /\* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */
+ /*   fclose(ficresfb); */
+ /*   printf("End of Computing Back forecasting \n"); */
+ /*   fprintf(ficlog,"End of Computing Back forecasting\n"); */
+ /* } */
  /************** Forecasting *****not tested NB*************/
  void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){
- Line 5484  void populforecast(char fileres[], doubl
+ Line 6550  void populforecast(char fileres[], doubl
    double calagedatem, agelim, kk1, kk2;
    double *popeffectif,*popcount;
    double ***p3mat,***tabpop,***tabpopprev;
-   double ***mobaverage;
+   /* double ***mobaverage; */
    char filerespop[FILENAMELENGTH];
    tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
- Line 5506  void populforecast(char fileres[], doubl
+ Line 6572  void populforecast(char fileres[], doubl
    if (cptcoveff==0) ncodemax[cptcoveff]=1;
-   if (mobilav!=0) {
+   /* if (mobilav!=0) { */
-     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+   /*   mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
-     if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
+   /*   if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ */
-       fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
+   /*     fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); */
-       printf(" Error in movingaverage mobilav=%d\n",mobilav);
+   /*     printf(" Error in movingaverage mobilav=%d\n",mobilav); */
-     }
+   /*   } */
-   }
+   /* } */
    stepsize=(int) (stepm+YEARM-1)/YEARM;
    if (stepm<=12) stepsize=1;
- Line 5521  void populforecast(char fileres[], doubl
+ Line 6587  void populforecast(char fileres[], doubl
    hstepm=1;
    hstepm=hstepm/stepm;
    if (popforecast==1) {
      if((ficpop=fopen(popfile,"r"))==NULL) {
        printf("Problem with population file : %s\n",popfile);exit(0);
- Line 5533  void populforecast(char fileres[], doubl
+ Line 6599  void populforecast(char fileres[], doubl
      i=1;
      while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF) i=i+1;
      imx=i;
      for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];
    }
    for(cptcov=1,k=0;cptcov<=i2;cptcov++){
-    for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
+     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
        k=k+1;
        fprintf(ficrespop,"\n#******");
        for(j=1;j<=cptcoveff;j++) {
- Line 5553  void populforecast(char fileres[], doubl
+ Line 6619  void populforecast(char fileres[], doubl
        for (cpt=0; cpt<=0;cpt++) {
          fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);
          for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){
            nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
            nhstepm = nhstepm/hstepm;
            p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
            oldm=oldms;savm=savms;
            hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
            for (h=0; h<=nhstepm; h++){
              if (h==(int) (calagedatem+YEARM*cpt)) {
                fprintf(ficrespop,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm);
- Line 5576  void populforecast(char fileres[], doubl
+ Line 6642  void populforecast(char fileres[], doubl
                }
                if (h==(int)(calagedatem+12*cpt)){
                  tabpop[(int)(agedeb)][j][cptcod]=kk1;
-                   /*fprintf(ficrespop," %.3f", kk1);
+                 /*fprintf(ficrespop," %.3f", kk1);
-                     if (popforecast==1) fprintf(ficrespop," [%.f]", kk1*popeffectif[(int)agedeb+1]);*/
+                   if (popforecast==1) fprintf(ficrespop," [%.f]", kk1*popeffectif[(int)agedeb+1]);*/
                }
              }
              for(i=1; i<=nlstate;i++){
                kk1=0.;
-                 for(j=1; j<=nlstate;j++){
+               for(j=1; j<=nlstate;j++){
-                   kk1= kk1+tabpop[(int)(agedeb)][j][cptcod];
+                 kk1= kk1+tabpop[(int)(agedeb)][j][cptcod];
-                 }
+               }
-                   tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedatem+12*cpt)*hstepm/YEARM*stepm-1)];
+               tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedatem+12*cpt)*hstepm/YEARM*stepm-1)];
              }
-             if (h==(int)(calagedatem+12*cpt)) for(j=1; j<=nlstate;j++)
+             if (h==(int)(calagedatem+12*cpt))
-               fprintf(ficrespop," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]);
+               for(j=1; j<=nlstate;j++)
+                 fprintf(ficrespop," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]);
            }
            free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
          }
        }
-   /******/
+       /******/
        for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) {
          fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);
          for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){
- Line 5621  void populforecast(char fileres[], doubl
+ Line 6688  void populforecast(char fileres[], doubl
            free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
          }
        }
-    }
+     }
    }
-   if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+   /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
    if (popforecast==1) {
      free_ivector(popage,0,AGESUP);
      free_vector(popeffectif,0,AGESUP);
- Line 5635  void populforecast(char fileres[], doubl
+ Line 6702  void populforecast(char fileres[], doubl
    free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
    fclose(ficrespop);
  } /* End of popforecast */
  int fileappend(FILE *fichier, char *optionfich)
  {
    if((fichier=fopen(optionfich,"a"))==NULL) {
- Line 5910  int readdata(char datafile[], int firsto
+ Line 6977  int readdata(char datafile[], int firsto
    if((fic=fopen(datafile,"r"))==NULL)    {
-     printf("Problem while opening datafile: %s\n", datafile);fflush(stdout);
+     printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);
-     fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);fflush(ficlog);return 1;
+     fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1;
    }
    i=1;
- Line 6288  int decodemodel ( char model[], int last
+ Line 7355  int decodemodel ( char model[], int last
  int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )
  {
    int i, m;
+   int firstone=0;
    for (i=1; i<=imx; i++) {
      for(m=2; (m<= maxwav); m++) {
        if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){
          anint[m][i]=9999;
-         s[m][i]=-1;
+         if (s[m][i] != -2) /* Keeping initial status of unknown vital status */
+           s[m][i]=-1;
        }
        if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){
          *nberr = *nberr + 1;
-         printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased (%d)\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr);
+         if(firstone == 0){
-         fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased (%d)\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr);
+           firstone=1;
+         printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results can be biased (%d) because status is a death state %d at wave %d. Wave dropped.\nOther similar cases in log file\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr,s[m][i],m);
+         }
+         fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results can be biased (%d) because status is a death state %d at wave %d. Wave dropped.\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr,s[m][i],m);
          s[m][i]=-1;
        }
        if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
- Line 6313  int calandcheckages(int imx, int maxwav,
+ Line 7385  int calandcheckages(int imx, int maxwav,
    for (i=1; i<=imx; i++)  {
      agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
      for(m=firstpass; (m<= lastpass); m++){
-       if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){
+       if(s[m][i] >0  || s[m][i]==-1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){ /* What if s[m][i]=-1 */
          if (s[m][i] >= nlstate+1) {
            if(agedc[i]>0){
              if((int)moisdc[i]!=99 && (int)andc[i]!=9999){
                agev[m][i]=agedc[i];
-           /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
+               /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
              }else {
                if ((int)andc[i]!=9999){
                  nbwarn++;
- Line 6328  int calandcheckages(int imx, int maxwav,
+ Line 7400  int calandcheckages(int imx, int maxwav,
                }
              }
            } /* agedc > 0 */
-         }
+         } /* end if */
          else if(s[m][i] !=9){ /* Standard case, age in fractional
                                   years but with the precision of a month */
            agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
- Line 6344  int calandcheckages(int imx, int maxwav,
+ Line 7416  int calandcheckages(int imx, int maxwav,
            }
            /*agev[m][i]=anint[m][i]-annais[i];*/
            /*     agev[m][i] = age[i]+2*m;*/
-         }
+         } /* en if 9*/
          else { /* =9 */
+           /* printf("Debug num[%d]=%ld s[%d][%d]=%d\n",i,num[i], m,i, s[m][i]); */
            agev[m][i]=1;
            s[m][i]=-1;
          }
        }
-       else /*= 0 Unknown */
+       else if(s[m][i]==0) /*= 0 Unknown */
          agev[m][i]=1;
-     }
+       else{
+         printf("Warning, num[%d]=%ld, s[%d][%d]=%d\n", i, num[i], m, i,s[m][i]);
+         fprintf(ficlog, "Warning, num[%d]=%ld, s[%d][%d]=%d\n", i, num[i], m, i,s[m][i]);
+         agev[m][i]=0;
+       }
+     } /* End for lastpass */
    }
    for (i=1; i<=imx; i++)  {
      for(m=firstpass; (m<=lastpass); m++){
        if (s[m][i] > (nlstate+ndeath)) {
- Line 6652  void syscompilerinfo(int logged)
+ Line 7730  void syscompilerinfo(int logged)
          return 0;
  }
+ int back_prevalence_limit(double *p, double **bprlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp, double dateprev1,double dateprev2, int firstpass, int lastpass, int mobilavproj){
+         /*--------------- Back Prevalence limit  (period or stable prevalence) --------------*/
+         /* Computes the back prevalence limit  for any combination      of covariate values
+    * at any age between ageminpar and agemaxpar
+          */
+   int i, j, k, i1 ;
+   /* double ftolpl = 1.e-10; */
+   double age, agebase, agelim;
+   double tot;
+   /* double ***mobaverage; */
+   /* double      **dnewm, **doldm, **dsavm;  /\* for use *\/ */
+   strcpy(fileresplb,"PLB_");
+   strcat(fileresplb,fileresu);
+   if((ficresplb=fopen(fileresplb,"w"))==NULL) {
+     printf("Problem with period (stable) back prevalence resultfile: %s\n", fileresplb);return 1;
+     fprintf(ficlog,"Problem with period (stable) back prevalence resultfile: %s\n", fileresplb);return 1;
+   }
+   printf("Computing period (stable) back prevalence: result on file '%s' \n", fileresplb);
+   fprintf(ficlog,"Computing period (stable) back prevalence: result on file '%s' \n", fileresplb);
+   pstamp(ficresplb);
+   fprintf(ficresplb,"# Period (stable) back prevalence. Precision given by ftolpl=%g \n", ftolpl);
+   fprintf(ficresplb,"#Age ");
+   for(i=1; i<=nlstate;i++) fprintf(ficresplb,"%d-%d ",i,i);
+   fprintf(ficresplb,"\n");
+   /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
+   agebase=ageminpar;
+   agelim=agemaxpar;
+   i1=pow(2,cptcoveff);
+   if (cptcovn < 1){i1=1;}
+   for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+     /* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
+     //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+     k=k+1;
+     /* to clean */
+     //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+     fprintf(ficresplb,"#******");
+     printf("#******");
+     fprintf(ficlog,"#******");
+     for(j=1;j<=cptcoveff;j++) {
+       fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+       printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+       fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+     }
+     fprintf(ficresplb,"******\n");
+     printf("******\n");
+     fprintf(ficlog,"******\n");
+     fprintf(ficresplb,"#Age ");
+     for(j=1;j<=cptcoveff;j++) {
+       fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+     }
+     for(i=1; i<=nlstate;i++) fprintf(ficresplb,"  %d-%d   ",i,i);
+     fprintf(ficresplb,"Total Years_to_converge\n");
+     for (age=agebase; age<=agelim; age++){
+       /* for (age=agebase; age<=agebase; age++){ */
+       if(mobilavproj > 0){
+         /* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */
+         /* bprevalim(bprlim, mobaverage, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
+         bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k);
+       }else if (mobilavproj == 0){
+         printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+         fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+         exit(1);
+       }else{
+         /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
+         bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);
+       }
+       fprintf(ficresplb,"%.0f ",age );
+       for(j=1;j<=cptcoveff;j++)
+         fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+       tot=0.;
+       for(i=1; i<=nlstate;i++){
+         tot +=  bprlim[i][i];
+         fprintf(ficresplb," %.5f", bprlim[i][i]);
+       }
+       fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
+     } /* Age */
+     /* was end of cptcod */
+   } /* cptcov */
+   /* hBijx(p, bage, fage); */
+   /* fclose(ficrespijb); */
+   return 0;
+ }
  int hPijx(double *p, int bage, int fage){
      /*------------- h Pij x at various ages ------------*/
- Line 6678  int hPijx(double *p, int bage, int fage)
+ Line 7852  int hPijx(double *p, int bage, int fage)
      agelim=AGESUP;
      hstepm=stepsize*YEARM; /* Every year of age */
      hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
      /* hstepm=1;   aff par mois*/
      pstamp(ficrespij);
      fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
      i1= pow(2,cptcoveff);
-    /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
+                 /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
-    /*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
+                 /*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
-    /*   k=k+1;  */
+                 /*      k=k+1;  */
      for (k=1; k <= (int) pow(2,cptcoveff); k++){
        fprintf(ficrespij,"\n#****** ");
        for(j=1;j<=cptcoveff;j++)
- Line 6719  int hPijx(double *p, int bage, int fage)
+ Line 7893  int hPijx(double *p, int bage, int fage)
        }
        /*}*/
      }
-         return 0;
+     return 0;
  }
+  int hBijx(double *p, int bage, int fage, double ***prevacurrent){
+     /*------------- h Bij x at various ages ------------*/
+   int stepsize;
+   /* int agelim; */
+         int ageminl;
+   int hstepm;
+   int nhstepm;
+   int h, i, i1, j, k;
+   double agedeb;
+   double ***p3mat;
+   strcpy(filerespijb,"PIJB_");  strcat(filerespijb,fileresu);
+   if((ficrespijb=fopen(filerespijb,"w"))==NULL) {
+     printf("Problem with Pij back resultfile: %s\n", filerespijb); return 1;
+     fprintf(ficlog,"Problem with Pij back resultfile: %s\n", filerespijb); return 1;
+   }
+   printf("Computing pij back: result on file '%s' \n", filerespijb);
+   fprintf(ficlog,"Computing pij back: result on file '%s' \n", filerespijb);
+   stepsize=(int) (stepm+YEARM-1)/YEARM;
+   /*if (stepm<=24) stepsize=2;*/
+   /* agelim=AGESUP; */
+   ageminl=30;
+   hstepm=stepsize*YEARM; /* Every year of age */
+   hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
+   /* hstepm=1;   aff par mois*/
+   pstamp(ficrespijb);
+   fprintf(ficrespijb,"#****** h Pij x Back Probability to be in state i at age x-h being in j at x ");
+   i1= pow(2,cptcoveff);
+   /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
+   /*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
+   /*    k=k+1;  */
+   for (k=1; k <= (int) pow(2,cptcoveff); k++){
+     fprintf(ficrespijb,"\n#****** ");
+     for(j=1;j<=cptcoveff;j++)
+       fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+     fprintf(ficrespijb,"******\n");
+     /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */
+     for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */
+       /* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /\* Typically 20 years = 20*12/6=40 *\/ */
+       nhstepm=(int) rint((agedeb-ageminl)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
+       nhstepm = nhstepm/hstepm; /* Typically 40/4=10, because estepm=24 stepm=6 => hstepm=24/6=4 */
+       /*          nhstepm=nhstepm*YEARM; aff par mois*/
+       p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+       /* oldm=oldms;savm=savms; */
+       /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);   */
+       hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k);
+       /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
+       fprintf(ficrespijb,"# Cov Agex agex-h hpijx with i,j=");
+       for(i=1; i<=nlstate;i++)
+         for(j=1; j<=nlstate+ndeath;j++)
+           fprintf(ficrespijb," %1d-%1d",i,j);
+       fprintf(ficrespijb,"\n");
+       for (h=0; h<=nhstepm; h++){
+         /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
+         fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb - h*hstepm/YEARM*stepm );
+         /* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */
+         for(i=1; i<=nlstate;i++)
+           for(j=1; j<=nlstate+ndeath;j++)
+             fprintf(ficrespijb," %.5f", p3mat[i][j][h]);
+         fprintf(ficrespijb,"\n");
+       }
+       free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+       fprintf(ficrespijb,"\n");
+     }
+     /*}*/
+   }
+   return 0;
+  } /*  hBijx */
  /***********************************************/
- Line 6753  int main(int argc, char *argv[])
+ Line 8004  int main(int argc, char *argv[])
    double agedeb=0.;
    double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW;
+         double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */
    double fret;
    double dum=0.; /* Dummy variable */
    double ***p3mat;
-   double ***mobaverage;
+   /* double ***mobaverage; */
    char line[MAXLINE];
    char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
- Line 6773  int main(int argc, char *argv[])
+ Line 8025  int main(int argc, char *argv[])
    int *tab;
    int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
+   int backcast=0;
    int mobilav=0,popforecast=0;
    int hstepm=0, nhstepm=0;
    int agemortsup;
- Line 6783  int main(int argc, char *argv[])
+ Line 8036  int main(int argc, char *argv[])
    double bage=0, fage=110., age, agelim=0., agebase=0.;
    double ftolpl=FTOL;
    double **prlim;
+   double **bprlim;
    double ***param; /* Matrix of parameters */
    double  *p;
    double **matcov; /* Matrix of covariance */
- Line 6794  int main(int argc, char *argv[])
+ Line 8048  int main(int argc, char *argv[])
    double *epj, vepp;
    double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
+   double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000;
    double **ximort;
    char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";
    int *dcwave;
- Line 7230  This is probably because your covariance
+ Line 8486  This is probably because your covariance
  Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model);
          exit(1);
        }else
-       if(mle==1)
+         if(mle==1)
-         printf("%1d%1d%1d",i1,j1,jk);
+           printf("%1d%1d%1d",i1,j1,jk);
        fprintf(ficlog,"%1d%1d%1d",i1,j1,jk);
        fprintf(ficparo,"%1d%1d%1d",i1,j1,jk);
        for(j=1; j <=i; j++){
- Line 7271  Please run with mle=-1 to get a correct
+ Line 8527  Please run with mle=-1 to get a correct
      }
      fprintf(ficres,"#%s\n",version);
    }    /* End of mle != -3 */
    /*  Main data
     */
    n= lastobs;
- Line 7353  Please run with mle=-1 to get a correct
+ Line 8609  Please run with mle=-1 to get a correct
    free_vector(annais,1,n);
    /* free_matrix(mint,1,maxwav,1,n);
       free_matrix(anint,1,maxwav,1,n);*/
-   free_vector(moisdc,1,n);
+   /* free_vector(moisdc,1,n); */
-   free_vector(andc,1,n);
+   /* free_vector(andc,1,n); */
    /* */
    wav=ivector(1,imx);
-   dh=imatrix(1,lastpass-firstpass+1,1,imx);
+   /* dh=imatrix(1,lastpass-firstpass+1,1,imx); */
-   bh=imatrix(1,lastpass-firstpass+1,1,imx);
+   /* bh=imatrix(1,lastpass-firstpass+1,1,imx); */
-   mw=imatrix(1,lastpass-firstpass+1,1,imx);
+   /* mw=imatrix(1,lastpass-firstpass+1,1,imx); */
+   dh=imatrix(1,lastpass-firstpass+2,1,imx); /* We are adding a wave if status is unknown at last wave but death occurs after last wave.*/
+   bh=imatrix(1,lastpass-firstpass+2,1,imx);
+   mw=imatrix(1,lastpass-firstpass+2,1,imx);
    /* Concatenates waves */
+   /* Concatenates waves: wav[i] is the number of effective (useful waves) of individual i.
+      Death is a valid wave (if date is known).
+      mw[mi][i] is the number of (mi=1 to wav[i]) effective wave out of mi of individual i
+      dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
+      and mw[mi+1][i]. dh depends on stepm.
+   */
    concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
    /* */
+   free_vector(moisdc,1,n);
+   free_vector(andc,1,n);
    /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */
    nbcode=imatrix(0,NCOVMAX,0,NCOVMAX);
- Line 7373  Please run with mle=-1 to get a correct
+ Line 8642  Please run with mle=-1 to get a correct
    Ndum =ivector(-1,NCOVMAX);
    if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */
      tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */
-   /* Nbcode gives the value of the lth modality of jth covariate, in
+   /* Nbcode gives the value of the lth modality (currently 1 to 2) of jth covariate, in
       V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/
-   /* 1 to ncodemax[j] is the maximum value of this jth covariate */
+   /* 1 to ncodemax[j] which is the maximum value of this jth covariate */
    /*  codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
    /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/
    /* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/
+   /* nbcode[Tvaraff[j]][codtabm(h,j)]) : if there are only 2 modalities for a covariate j,
+    * codtabm(h,j) gives its value classified at position h and nbcode gives how it is coded
+    * (currently 0 or 1) in the data.
+    * In a loop on h=1 to 2**k, and a loop on j (=1 to k), we get the value of
+    * corresponding modality (h,j).
+    */
    h=0;
- Line 7389  Please run with mle=-1 to get a correct
+ Line 8665  Please run with mle=-1 to get a correct
    m=pow(2,cptcoveff);
            /**< codtab(h,k)  k   = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1
-            * For k=4 covariates, h goes from 1 to 2**k
+            * For k=4 covariates, h goes from 1 to m=2**k
-            * codtabm(h,k)=  1 & (h-1) >> (k-1) ;
+            * codtabm(h,k)=  (1 & (h-1) >> (k-1)) + 1;
+            * #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
             *     h\k   1     2     3     4
             *______________________________
             *     1 i=1 1 i=1 1 i=1 1 i=1 1
- Line 7410  Please run with mle=-1 to get a correct
+ Line 8687  Please run with mle=-1 to get a correct
             *    15 i=8 1     2     2     2
             *    16     2     2     2     2
             */
+   /* How to do the opposite? From combination h (=1 to 2**k) how to get the value on the covariates? */
+      /* from h=5 and m, we get then number of covariates k=log(m)/log(2)=4
+      * and the value of each covariate?
+      * V1=1, V2=1, V3=2, V4=1 ?
+      * h-1=4 and 4 is 0100 or reverse 0010, and +1 is 1121 ok.
+      * h=6, 6-1=5, 5 is 0101, 1010, 2121, V1=2nd, V2=1st, V3=2nd, V4=1st.
+      * In order to get the real value in the data, we use nbcode
+      * nbcode[Tvar[3][2nd]]=1 and nbcode[Tvar[4][1]]=0
+      * We are keeping this crazy system in order to be able (in the future?)
+      * to have more than 2 values (0 or 1) for a covariate.
+      * #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
+      * h=6, k=2? h-1=5=0101, reverse 1010, +1=2121, k=2nd position: value is 1: codtabm(6,2)=1
+      *              bbbbbbbb
+      *              76543210
+      *   h-1        00000101 (6-1=5)
+      *(h-1)>>(k-1)= 00000001 >> (2-1) = 1 right shift
+      *           &
+      *     1        00000001 (1)
+      *              00000001        = 1 & ((h-1) >> (k-1))
+      *          +1= 00000010 =2
+      *
+      * h=14, k=3 => h'=h-1=13, k'=k-1=2
+      *          h'      1101 =2^3+2^2+0x2^1+2^0
+      *    >>k'            11
+      *          &   00000001
+      *            = 00000001
+      *      +1    = 00000010=2    =  codtabm(14,3)
+      * Reverse h=6 and m=16?
+      * cptcoveff=log(16)/log(2)=4 covariate: 6-1=5=0101 reversed=1010 +1=2121 =>V1=2, V2=1, V3=2, V4=1.
+      * for (j=1 to cptcoveff) Vj=decodtabm(j,h,cptcoveff)
+      * decodtabm(h,j,cptcoveff)= (((h-1) >> (j-1)) & 1) +1
+      * decodtabm(h,j,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (j-1)) & 1) +1 : -1)
+      * V3=decodtabm(14,3,2**4)=2
+      *          h'=13   1101 =2^3+2^2+0x2^1+2^0
+      *(h-1) >> (j-1)    0011 =13 >> 2
+      *          &1 000000001
+      *           = 000000001
+      *         +1= 000000010 =2
+      *                  2211
+      *                  V1=1+1, V2=0+1, V3=1+1, V4=1+1
+      *                  V3=2
+      */
    /* /\* for(h=1; h <=100 ;h++){  *\/ */
    /*   /\* printf("h=%2d ", h); *\/ */
    /*    /\* for(k=1; k <=10; k++){ *\/ */
- Line 7488  Title=%s <br>Datafile=%s Firstpass=%d La
+ Line 8808  Title=%s <br>Datafile=%s Firstpass=%d La
            optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
    }
-   fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C)  2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015</a></font><br>  \
+   fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C)  2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-2013-2016-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br>  \
  <hr size=\"2\" color=\"#EC5E5E\"> \n\
  <font size=\"2\">IMaCh-%s <br> %s</font> \
  <hr size=\"2\" color=\"#EC5E5E\"> \n\
- Line 7518  Title=%s <br>Datafile=%s Firstpass=%d La
+ Line 8838  Title=%s <br>Datafile=%s Firstpass=%d La
    /* Calculates basic frequencies. Computes observed prevalence at single age
       and prints on file fileres'p'. */
-   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart);
+   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\
+               firstpass, lastpass,  stepm,  weightopt, model);
    fprintf(fichtm,"\n");
    fprintf(fichtm,"<br>Total number of observations=%d <br>\n\
- Line 7530  Interval (in months) between two waves:
+ Line 8851  Interval (in months) between two waves:
      newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
      savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
      oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
    /* For Powell, parameters are in a vector p[] starting at p[1]
       so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
    p=param[1][1]; /* *(*(*(param +1)+1)+0) */
- Line 7549  Interval (in months) between two waves:
+ Line 8869  Interval (in months) between two waves:
      for (i=1; i<=imx; i++){
        dcwave[i]=-1;
        for (m=firstpass; m<=lastpass; m++)
-         if (s[m][i]>nlstate) {
+                                 if (s[m][i]>nlstate) {
-           dcwave[i]=m;
+                                         dcwave[i]=m;
-           /*    printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/
+                                         /*      printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/
-           break;
+                                         break;
-         }
+                                 }
      }
      for (i=1; i<=imx; i++) {
        if (wav[i]>0){
-         ageexmed[i]=agev[mw[1][i]][i];
+                                 ageexmed[i]=agev[mw[1][i]][i];
-         j=wav[i];
+                                 j=wav[i];
-         agecens[i]=1.;
+                                 agecens[i]=1.;
-         if (ageexmed[i]> 1 && wav[i] > 0){
+                                 if (ageexmed[i]> 1 && wav[i] > 0){
-           agecens[i]=agev[mw[j][i]][i];
+                                         agecens[i]=agev[mw[j][i]][i];
-           cens[i]= 1;
+                                         cens[i]= 1;
-         }else if (ageexmed[i]< 1)
+                                 }else if (ageexmed[i]< 1)
-           cens[i]= -1;
+                                         cens[i]= -1;
-         if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass)
+                                 if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass)
-           cens[i]=0 ;
+                                         cens[i]=0 ;
        }
        else cens[i]=-1;
      }
      for (i=1;i<=NDIM;i++) {
        for (j=1;j<=NDIM;j++)
-         ximort[i][j]=(i == j ? 1.0 : 0.0);
+                                 ximort[i][j]=(i == j ? 1.0 : 0.0);
      }
      /*p[1]=0.0268; p[NDIM]=0.083;*/
- Line 7981  Please run with mle=-1 to get a correct
+ Line 9301  Please run with mle=-1 to get a correct
      }
      fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
-     fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
+     fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
-     fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
+     fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d, ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
      /* Other stuffs, more or less useful */
      while((c=getc(ficpar))=='#' && c!= EOF){
- Line 8031  Please run with mle=-1 to get a correct
+ Line 9351  Please run with mle=-1 to get a correct
      fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
      /* day and month of proj2 are not used but only year anproj2.*/
+     while((c=getc(ficpar))=='#' && c!= EOF){
+       ungetc(c,ficpar);
+       fgets(line, MAXLINE, ficpar);
+       fputs(line,stdout);
+       fputs(line,ficparo);
+     }
+     ungetc(c,ficpar);
+     fscanf(ficpar,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
+     fscanf(ficparo,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
+     fscanf(ficlog,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
+     fscanf(ficres,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
+     /* day and month of proj2 are not used but only year anproj2.*/
       /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */
- Line 8045  Please run with mle=-1 to get a correct
+ Line 9378  Please run with mle=-1 to get a correct
  This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
  Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
      }else
-       printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
+       printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p);
      printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\
-                  model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\
+                  model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,backcast, estepm, \
-                  jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
+                  jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2);
     /*------------ free_vector  -------------*/
     /*  chdir(path); */
-     free_ivector(wav,1,imx);
+     /* free_ivector(wav,1,imx); */  /* Moved after last prevalence call */
-     free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
+     /* free_imatrix(dh,1,lastpass-firstpass+2,1,imx); */
-     free_imatrix(bh,1,lastpass-firstpass+1,1,imx);
+     /* free_imatrix(bh,1,lastpass-firstpass+2,1,imx); */
-     free_imatrix(mw,1,lastpass-firstpass+1,1,imx);
+     /* free_imatrix(mw,1,lastpass-firstpass+2,1,imx);    */
      free_lvector(num,1,n);
      free_vector(agedc,1,n);
      /*free_matrix(covar,0,NCOVMAX,1,n);*/
- Line 8075  Please run with mle=-1 to get a correct
+ Line 9408  Please run with mle=-1 to get a correct
      prevalence_limit(p, prlim,  ageminpar, agemaxpar, ftolpl, &ncvyear);
      fclose(ficrespl);
- #ifdef FREEEXIT2
- #include "freeexit2.h"
- #endif
      /*------------- h Pij x at various ages ------------*/
      /*#include "hpijx.h"*/
      hPijx(p, bage, fage);
      fclose(ficrespij);
+                 ncovcombmax=  pow(2,cptcoveff);
    /*-------------- Variance of one-step probabilities---*/
      k=1;
      varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
+                 /* Prevalence for each covariates in probs[age][status][cov] */
-     probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
+     probs= ma3x(1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
      for(i=1;i<=AGESUP;i++)
-       for(j=1;j<=NCOVMAX;j++)
+       for(j=1;j<=nlstate;j++)
-         for(k=1;k<=NCOVMAX;k++)
+                                 for(k=1;k<=ncovcombmax;k++)
-           probs[i][j][k]=0.;
+                                         probs[i][j][k]=0.;
+                 prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+                 if (mobilav!=0 ||mobilavproj !=0 ) {
+                         mobaverage= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
+                         if (mobilav!=0) {
+                                 if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilav)!=0){
+                                         fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
+                                         printf(" Error in movingaverage mobilav=%d\n",mobilav);
+                                 }
+                         }
+                         /* /\* Prevalence for each covariates in probs[age][status][cov] *\/ */
+                         /* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
+                         else if (mobilavproj !=0) {
+                                 if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){
+                                         fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj);
+                                         printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj);
+                                 }
+                         }
+                 }/* end if moving average */
      /*---------- Forecasting ------------------*/
      /*if((stepm == 1) && (strcmp(model,".")==0)){*/
      if(prevfcast==1){
        /*    if(stepm ==1){*/
        prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
-       /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/
-       /*      }  */
-       /*      else{ */
-       /*        erreur=108; */
-       /*        printf("Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
-       /*        fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
-       /*      } */
      }
+     if(backcast==1){
+                 ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
+                 ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
+                 ddsavms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
+                         /*--------------- Back Prevalence limit  (period or stable prevalence) --------------*/
+                         /*#include "prevlim.h"*/  /* Use ficresplb, ficlog */
+                         bprlim=matrix(1,nlstate,1,nlstate);
+                         back_prevalence_limit(p, bprlim,  ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj);
+                         fclose(ficresplb);
+                         hBijx(p, bage, fage, mobaverage);
+                         fclose(ficrespijb);
+                         free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
+       /* prevbackforecast(fileresu, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff); */
+                         free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath);
+                         free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath);
+                         free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath);
+                 }
+   /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+   /* if (mobilavproj!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+     /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/
+     /*      }  */
+     /*      else{ */
+     /*        erreur=108; */
+     /*        printf("Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
+     /*        fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
+     /*      } */
      /* ------ Other prevalence ratios------------ */
      /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */
-     prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+     /* prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
      /*  printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d,  mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\
          ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);
      */
+     free_ivector(wav,1,imx);
-     if (mobilav!=0) {
+     free_imatrix(dh,1,lastpass-firstpass+2,1,imx);
-       mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+     free_imatrix(bh,1,lastpass-firstpass+2,1,imx);
-       if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
+     free_imatrix(mw,1,lastpass-firstpass+2,1,imx);
-         fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
-         printf(" Error in movingaverage mobilav=%d\n",mobilav);
-       }
+                 /*   if (mobilav!=0) { */
-     }
+                 /*     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+                 /*     if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ */
+                 /* fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); */
+                 /* printf(" Error in movingaverage mobilav=%d\n",mobilav); */
+                 /*     } */
+                 /*   } */
      /*---------- Health expectancies, no variances ------------*/
      strcpy(filerese,"E_");
      strcat(filerese,fileresu);
      if((ficreseij=fopen(filerese,"w"))==NULL) {
- Line 8139  Please run with mle=-1 to get a correct
+ Line 9516  Please run with mle=-1 to get a correct
      fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog);
      /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
        for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
      for (k=1; k <= (int) pow(2,cptcoveff); k++){
-         fprintf(ficreseij,"\n#****** ");
+                         fprintf(ficreseij,"\n#****** ");
-         for(j=1;j<=cptcoveff;j++) {
+                         for(j=1;j<=cptcoveff;j++) {
-           fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+                                 fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-         }
+                         }
-         fprintf(ficreseij,"******\n");
+                         fprintf(ficreseij,"******\n");
-         eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
+                         eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
-         oldm=oldms;savm=savms;
+                         oldm=oldms;savm=savms;
-         evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);
+                         evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);
-         free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
+                         free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
        /*}*/
      }
      fclose(ficreseij);
      printf("done evsij\n");fflush(stdout);
      fprintf(ficlog,"done evsij\n");fflush(ficlog);
      /*---------- Health expectancies and variances ------------*/
      strcpy(filerest,"T_");
      strcat(filerest,fileresu);
      if((ficrest=fopen(filerest,"w"))==NULL) {
- Line 8169  Please run with mle=-1 to get a correct
+ Line 9546  Please run with mle=-1 to get a correct
      }
      printf("Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(stdout);
      fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(ficlog);
      strcpy(fileresstde,"STDE_");
      strcat(fileresstde,fileresu);
- Line 8204  Please run with mle=-1 to get a correct
+ Line 9581  Please run with mle=-1 to get a correct
      for (k=1; k <= (int) pow(2,cptcoveff); k++){
        fprintf(ficrest,"\n#****** ");
        for(j=1;j<=cptcoveff;j++)
-         fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+                                 fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
        fprintf(ficrest,"******\n");
        fprintf(ficresstdeij,"\n#****** ");
        fprintf(ficrescveij,"\n#****** ");
        for(j=1;j<=cptcoveff;j++) {
-         fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+                                 fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-         fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+                                 fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
        }
        fprintf(ficresstdeij,"******\n");
        fprintf(ficrescveij,"******\n");
        fprintf(ficresvij,"\n#****** ");
        for(j=1;j<=cptcoveff;j++)
-         fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+                                 fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
        fprintf(ficresvij,"******\n");
        eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
- Line 8238  Please run with mle=-1 to get a correct
+ Line 9615  Please run with mle=-1 to get a correct
        for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
-         oldm=oldms;savm=savms; /* ZZ Segmentation fault */
+                                 oldm=oldms;savm=savms; /* ZZ Segmentation fault */
-         cptcod= 0; /* To be deleted */
+                                 cptcod= 0; /* To be deleted */
-         printf("varevsij %d \n",vpopbased);
+                                 printf("varevsij %d \n",vpopbased);
-         fprintf(ficlog, "varevsij %d \n",vpopbased);
+                                 fprintf(ficlog, "varevsij %d \n",vpopbased);
-         varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
+                                 varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
-         fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");
+                                 fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");
-         if(vpopbased==1)
+                                 if(vpopbased==1)
-           fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
+                                         fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
-         else
+                                 else
-           fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
+                                         fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
-         fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
+                                 fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
-         for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
+                                 for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
-         fprintf(ficrest,"\n");
+                                 fprintf(ficrest,"\n");
-         /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
+                                 /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
-         epj=vector(1,nlstate+1);
+                                 epj=vector(1,nlstate+1);
-         printf("Computing age specific period (stable) prevalences in each health state \n");
+                                 printf("Computing age specific period (stable) prevalences in each health state \n");
-         fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
+                                 fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
-         for(age=bage; age <=fage ;age++){
+                                 for(age=bage; age <=fage ;age++){
-           prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */
+                                         prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */
-           if (vpopbased==1) {
+                                         if (vpopbased==1) {
-             if(mobilav ==0){
+                                                 if(mobilav ==0){
-               for(i=1; i<=nlstate;i++)
+                                                         for(i=1; i<=nlstate;i++)
-                 prlim[i][i]=probs[(int)age][i][k];
+                                                                 prlim[i][i]=probs[(int)age][i][k];
-             }else{ /* mobilav */
+                                                 }else{ /* mobilav */
-               for(i=1; i<=nlstate;i++)
+                                                         for(i=1; i<=nlstate;i++)
-                 prlim[i][i]=mobaverage[(int)age][i][k];
+                                                                 prlim[i][i]=mobaverage[(int)age][i][k];
-             }
+                                                 }
-           }
+                                         }
-           fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
+                                         fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
-           /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
+                                         /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
-           /* printf(" age %4.0f ",age); */
+                                         /* printf(" age %4.0f ",age); */
-           for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
+                                         for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
-             for(i=1, epj[j]=0.;i <=nlstate;i++) {
+                                                 for(i=1, epj[j]=0.;i <=nlstate;i++) {
-               epj[j] += prlim[i][i]*eij[i][j][(int)age];
+                                                         epj[j] += prlim[i][i]*eij[i][j][(int)age];
-               /*ZZZ  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+                                                         /*ZZZ  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
-               /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
+                                                         /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
-             }
+                                                 }
-             epj[nlstate+1] +=epj[j];
+                                                 epj[nlstate+1] +=epj[j];
-           }
+                                         }
-           /* printf(" age %4.0f \n",age); */
+                                         /* printf(" age %4.0f \n",age); */
-           for(i=1, vepp=0.;i <=nlstate;i++)
+                                         for(i=1, vepp=0.;i <=nlstate;i++)
-             for(j=1;j <=nlstate;j++)
+                                                 for(j=1;j <=nlstate;j++)
-               vepp += vareij[i][j][(int)age];
+                                                         vepp += vareij[i][j][(int)age];
-           fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
+                                         fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
-           for(j=1;j <=nlstate;j++){
+                                         for(j=1;j <=nlstate;j++){
-             fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
+                                                 fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
-           }
+                                         }
-           fprintf(ficrest,"\n");
+                                         fprintf(ficrest,"\n");
-         }
+                                 }
        } /* End vpopbased */
        free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
        free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
- Line 8329  Please run with mle=-1 to get a correct
+ Line 9706  Please run with mle=-1 to get a correct
      for (k=1; k <= (int) pow(2,cptcoveff); k++){
          fprintf(ficresvpl,"\n#****** ");
-         for(j=1;j<=cptcoveff;j++)
+                         for(j=1;j<=cptcoveff;j++)
-           fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+                                 fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-         fprintf(ficresvpl,"******\n");
+                         fprintf(ficresvpl,"******\n");
-         varpl=matrix(1,nlstate,(int) bage, (int) fage);
+                         varpl=matrix(1,nlstate,(int) bage, (int) fage);
-         oldm=oldms;savm=savms;
+                         oldm=oldms;savm=savms;
-         varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart);
+                         varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart);
-         free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
+                         free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
        /*}*/
      }
      fclose(ficresvpl);
      printf("done variance-covariance of period prevalence\n");fflush(stdout);
      fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);
      /*---------- End : free ----------------*/
-     if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+     if (mobilav!=0 ||mobilavproj !=0) free_ma3x(mobaverage,1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */
-     free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
+     free_ma3x(probs,1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
    }  /* mle==-3 arrives here for freeing */
   /* endfree:*/
      free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */
- Line 8376  Please run with mle=-1 to get a correct
+ Line 9753  Please run with mle=-1 to get a correct
    if((nberr >0) || (nbwarn>0)){
-     printf("End of Imach with %d errors and/or %d warnings\n",nberr,nbwarn);
+     printf("End of Imach with %d errors and/or %d warnings. Please look at the log file for details.\n",nberr,nbwarn);
-     fprintf(ficlog,"End of Imach with %d errors and/or warnings %d\n",nberr,nbwarn);
+     fprintf(ficlog,"End of Imach with %d errors and/or warnings %d. Please look at the log file for details.\n",nberr,nbwarn);
    }else{
      printf("End of Imach\n");
      fprintf(ficlog,"End of Imach\n");

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