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.215 and 1.219

-version 1.215, 2015/12/16 08:52:24
+version 1.219, 2016/02/15 00:48:12
  Line 1
  /* $Id$
    $State$
    $Log$
+   Revision 1.219  2016/02/15 00:48:12  brouard
+   *** empty log message ***
+   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
- Line 640
+ Line 657
    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 778  typedef struct {
+ Line 812  typedef struct {
  #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 810  int cptcovage=0; /**< Number of covariat
+ Line 846  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 832  double jmean=1; /* Mean space between 2
+ Line 869  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, *ficresphtm, *ficresphtmfr, *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 856  char fileresv[FILENAMELENGTH];
+ Line 895  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 939  int *ncodemaxwundef;  /* ncodemax[j]= Nu
+ Line 978  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, ***mobaverages; /* New global variable */
  double *ageexmed,*agecens;
  double dateintmean=0;
- Line 2022  void powell(double p[], double **xi, int
+ Line 2062  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 2045  double **prevalim(double **prlim, int nl
+ Line 2085  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 2054  double **prevalim(double **prlim, int nl
+ Line 2094  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 2072  double **prevalim(double **prlim, int nl
+ Line 2113  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 2087  double **prevalim(double **prlim, int nl
+ Line 2129  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 2133  Earliest age to start was %d-%d=%d, ncvl
+ Line 2176  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 2154  double **pmij(double **ps, double *cov,
+ Line 2328  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 2239  double **matprod2(double **out, double *
+ Line 2573  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 2272  double ***hpxij(double ***po, int nhstep
+ Line 2606  double ***hpxij(double ***po, int nhstep
        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];
+                                 /*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;
+ }
+ /************* 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];
+                                 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;
  }
  #ifdef NLOPT
    double  myfunc(unsigned n, const double *p1, double *grad, void *pd){
    double fret;
- Line 2446  double func( double *x)
+ Line 2885  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 2578  double func( double *x)
+ Line 3014  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 2685  double funcone( double *x)
+ Line 3125  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 3070  double hessij( double x[], double **hess
+ Line 3514  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 3242  void  freqsummary(char fileres[], int ia
+ Line 3686  void  freqsummary(char fileres[], int ia
    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);*/
- Line 3282  Title=%s <br>Datafile=%s Firstpass=%d La
+ Line 3727  Title=%s <br>Datafile=%s Firstpass=%d La
    }
    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,iagemax+3);
+   freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
    j1=0;
    j=cptcoveff;
- Line 3491  Title=%s <br>Datafile=%s Firstpass=%d La
+ Line 3936  Title=%s <br>Datafile=%s Firstpass=%d La
    fclose(ficresp);
    fclose(ficresphtm);
    fclose(ficresphtmfr);
-   free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3);
+   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 3519  void prevalence(double ***probs, double
+ Line 3964  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 3527  void prevalence(double ***probs, double
+ Line 3972  void prevalence(double ***probs, double
    if (cptcovn<1) {j=1;ncodemax[1]=1;}
    first=1;
-   for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){
+   for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */
      for (i=1; i<=nlstate; i++)
-       for(iage=iagemin; iage <= iagemax+3; iage++)
+       for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++)
-         prop[i][iage]=0.0;
+                                 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++)
+                                 for (z1=1; z1<=cptcoveff; z1++) /* For each covariate, look at the value for individual i and checks if it is equal to the corresponding value of this covariate according to current combination j1*/
-           if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)])
+                                         if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)])
-             bool=0;
+                                                 bool=0;
        }
-       if (bool==1) {
+       if (bool==1) { /* For this combination of covariates values, this individual fits */
-         /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */
+                                 /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */
-         for(mi=1; mi<wav[i];mi++){
+                                 for(mi=1; mi<wav[i];mi++){
-           m=mw[mi][i];
+                                         m=mw[mi][i];
-           agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
+                                         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 *\/ */
+                                         /* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */
-           if(m >=firstpass && m <=lastpass){
+                                         if(m >=firstpass && m <=lastpass){
-             y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */
+                                                 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 ((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]==0) agev[m][i]=iagemax+1;
-               if(agev[m][i]==1) agev[m][i]=iagemax+2;
+                                                         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((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+3+AGEMARGE){
-               if (s[m][i]>0 && s[m][i]<=nlstate) {
+                                                                 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);
-                 /*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]]);*/
+                                                                 exit(1);
-                 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];
+                                                         if (s[m][i]>0 && s[m][i]<=nlstate) {
-               } /* end valid statuses */
+                                                                 /*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]]);*/
-             } /* end selection of dates */
+                                                                 prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */
-           } /* end selection of waves */
+                                                                 prop[s[m][i]][iagemax+3] += weight[i];
-         } /* end effective waves */
+                                                         } /* 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];
+                                 posprop += prop[jk][i];
        }
        for(jk=1; jk <=nlstate ; jk++){
-         if( i <=  iagemax){
+                                 if( i <=  iagemax){
-           if(posprop>=1.e-5){
+                                         if(posprop>=1.e-5){
-             probs[i][jk][j1]= prop[jk][i]/posprop;
+                                                 probs[i][jk][j1]= prop[jk][i]/posprop;
-           } else{
+                                         } else{
-             if(first==1){
+                                                 if(first==1){
-               first=0;
+                                                         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]);
+                                                         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 */
- Line 3584  void prevalence(double ***probs, double
+ Line 4032  void prevalence(double ***probs, double
    /*  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 3601  void  concatwav(int wav[], int **dh, int
+ Line 4049  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, firstwo;
+   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.;
- Line 3613  void  concatwav(int wav[], int **dh, int
+ Line 4063  void  concatwav(int wav[], int **dh, int
      mi=0;
      m=firstpass;
      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 */
+     /* 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 (andc[i] != 9999) { /* A death occured after lastpass */
+     }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++;
+       /* m++; */
-       mi++;
+       /* mi++; */
-       s[m][i]=nlstate+1;  /* We are setting the status to the last of non live state */
+       /* s[m][i]=nlstate+1;  /\* We are setting the status to the last of non live state *\/ */
-       mw[mi][i]=m;
+       /* mw[mi][i]=m; */
-       nbwarn++;
+       nberr++;
-       if(firstwo==0){
+       if ((int)anint[m][i]!= 9999) { /* date of last interview is known */
-         printf("Warning! Death for individual %ld line=%d  occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\nOthers in log file only\n",num[i],i,lastpass,nlstate+1, m);
+         if(firstwo==0){
-         fprintf(ficlog,"Warning! Death for individual %ld line=%d  occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\n",num[i],i,lastpass,nlstate+1, m);
+           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;
+           firstwo=1;
-       }
+         }
-       if(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 );
-         fprintf(ficlog,"Warning! Death for individual %ld line=%d  occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\n",num[i],i,lastpass,nlstate+1, 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;
- Line 3769  void tricode(int *Tvar, int **nbcode, in
+ Line 4245  void tricode(int *Tvar, int **nbcode, in
    for (j=1; j<=(cptcovs); j++) { /* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only */
      for (k=-1; k < maxncov; k++) Ndum[k]=0;
      for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the
-                                modality of this covariate Vj*/
+                                                                                                                                 modality of this covariate Vj*/
        ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
-                                     * If product of Vn*Vm, still boolean *:
+                                                                                                                                                 * If product of Vn*Vm, still boolean *:
-                                     * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
+                                                                                                                                                 * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
-                                     * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0   */
+                                                                                                                                                 * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0   */
        /* Finds for covariate j, n=Tvar[j] of Vn . ij is the
                                        modality of the nth covariate of individual i. */
        if (ij > modmaxcovj)
          modmaxcovj=ij;
        else if (ij < modmincovj)
-         modmincovj=ij;
+                                 modmincovj=ij;
        if ((ij < -1) && (ij > NCOVMAX)){
-         printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
+                                 printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
-         exit(1);
+                                 exit(1);
        }else
        Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/
        /*  If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */
- Line 3800  void tricode(int *Tvar, int **nbcode, in
+ Line 4276  void tricode(int *Tvar, int **nbcode, in
        printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]);
        fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]);
        if( Ndum[k] != 0 ){ /* Counts if nobody answered modality k ie empty modality, we skip it and reorder */
-         if( k != -1){
+                                 if( k != -1){
-           ncodemax[j]++;  /* ncodemax[j]= Number of modalities of the j th
+                                         ncodemax[j]++;  /* ncodemax[j]= Number of modalities of the j th
-                              covariate for which somebody answered excluding
+                                                                                                                  covariate for which somebody answered excluding
-                              undefined. Usually 2: 0 and 1. */
+                                                                                                                  undefined. Usually 2: 0 and 1. */
-         }
+                                 }
-         ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th
+                                 ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th
-                              covariate for which somebody answered including
+                                                                                                                                 covariate for which somebody answered including
-                              undefined. Usually 3: -1, 0 and 1. */
+                                                                                                                                 undefined. Usually 3: -1, 0 and 1. */
        }
        /* In fact  ncodemax[j]=2 (dichotom. variables only) but it could be more for
-          historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
+                                  historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
      } /* Ndum[-1] number of undefined modalities */
      /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
      /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7.
         If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125;
- Line 3829  void tricode(int *Tvar, int **nbcode, in
+ Line 4305  void tricode(int *Tvar, int **nbcode, in
      ij=0; /* ij is similar to i but can jump over null modalities */
      for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/
          if (Ndum[i] == 0) { /* If nobody responded to this modality k */
-           break;
+                                 break;
-         }
+                         }
          ij++;
          nbcode[Tvar[j]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
          cptcode = ij; /* New max modality for covar j */
- Line 3851  void tricode(int *Tvar, int **nbcode, in
+ Line 4327  void tricode(int *Tvar, int **nbcode, in
      /*   }  /\* end of loop on modality k *\/ */
    } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/
-  for (k=-1; k< maxncov; k++) Ndum[k]=0;
+         for (k=-1; k< maxncov; k++) Ndum[k]=0;
    for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */
-    /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
+                 /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
-    ij=Tvar[i]; /* Tvar might be -1 if status was unknown */
+                 ij=Tvar[i]; /* Tvar might be -1 if status was unknown */
-    Ndum[ij]++; /* Might be supersed V1 + V1*age */
+                 Ndum[ij]++; /* Might be supersed V1 + V1*age */
-  }
+         }
-  ij=0;
+         ij=0;
-  for (i=0; i<=  maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
+         for (i=0; i<=  maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
-    /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
+                 /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
-    if((Ndum[i]!=0) && (i<=ncovcol)){
+                 if((Ndum[i]!=0) && (i<=ncovcol)){
-      ij++;
+                         ij++;
-      /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
+                         /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
-      Tvaraff[ij]=i; /*For printing (unclear) */
+                         Tvaraff[ij]=i; /*For printing (unclear) */
-    }else{
+                 }else{
-        /* Tvaraff[ij]=0; */
+                         /* Tvaraff[ij]=0; */
-    }
+                 }
-  }
+         }
-  /* ij--; */
+         /* ij--; */
-  cptcoveff=ij; /*Number of total covariates*/
+         cptcoveff=ij; /*Number of total covariates*/
  }
- Line 4084  void cvevsij(double ***eij, double x[],
+ Line 4560  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 4187  void cvevsij(double ***eij, double x[],
+ Line 4663  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 4969  To be simple, these graphs help to under
+ Line 5445  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 prevfcast, int estepm ,          \
+                   int popforecast, int prevfcast, int backcast, int estepm , \
                    double jprev1, double mprev1,double anprev1, double dateprev1, \
                    double jprev2, double mprev2,double anprev2, double dateprev2){
    int jj1, k1, i1, cpt;
- Line 4987  void printinghtml(char fileresu[], char
+ Line 5463  void printinghtml(char fileresu[], char
   - 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,"\
+  - 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_"));
- Line 5041  divided by h: <sub>h</sub>P<sub>ij</sub>
+ Line 5523  divided by h: <sub>h</sub>P<sub>ij</sub>
       }
       /* 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++){
- Line 5122  See page 'Matrix of variance-covariance
+ Line 5611  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 5139  true period expectancies (those weighted
+ Line 5628  true period expectancies (those weighted
  }
  /******************* Gnuplot file **************/
-     void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, 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;
+         int ioffset; /* variable offset for columns */
  /*   if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
  /*     printf("Problem with file %s",optionfilegnuplot); */
  /*     fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */
- Line 5175  true period expectancies (those weighted
+ Line 5666  true period expectancies (those weighted
        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:($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:($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,",\\\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 5192  true period expectancies (those weighted
+ Line 5683  true period expectancies (those weighted
      for (k1=1; k1<= m ; k1 ++) { /* For each combination of covariate */
        /* 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 */
+       for (k=1; k<=cptcoveff; k++){    /* For each covariate k get corresponding value lv for combination k1 */
-         lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+                                 lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
-         /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+                                 /* 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(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 */
+                                 /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-         vlv= nbcode[Tvaraff[lv]][lv];
+                                 vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */
-         fprintf(ficgp," V%d=%d ",k,vlv);
+                         /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
+                                 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,"\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,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
-      fprintf(ficgp,"set xlabel \"Age\" \n\
+                         fprintf(ficgp,"set xlabel \"Age\" \n\
- set ylabel \"Probability\" \n\
+ set ylabel \"Probability\" \n   \
- set ter svg size 640, 480\n\
+ set ter svg size 640, 480\n     \
  plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
-      for (i=1; i<= nlstate ; i ++) {
+                         for (i=1; i<= nlstate ; i ++) {
-        if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+                                 if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
-        else        fprintf(ficgp," %%*lf (%%*lf)");
+                                 else        fprintf(ficgp," %%*lf (%%*lf)");
-      }
+                         }
-      fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+                         fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
-      for (i=1; i<= nlstate ; i ++) {
+                         for (i=1; i<= nlstate ; i ++) {
-        if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+                                 if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
-        else fprintf(ficgp," %%*lf (%%*lf)");
+                                 else fprintf(ficgp," %%*lf (%%*lf)");
-      }
+                         }
-      fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+                         fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
-      for (i=1; i<= nlstate ; i ++) {
+                         for (i=1; i<= nlstate ; i ++) {
-        if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+                                 if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
-        else 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));
-      fprintf(ficgp,"\nset out \n");
+                         if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */
+                                 /* 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,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */
+                                 kl=0;
+                                 for (k=1; k<=cptcoveff; k++){    /* For each combination of 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[k]][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){
+                                                         fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' with line ",kl+1, k,kl+1+1,nbcode[Tvaraff[k]][lv], \
++(cpt-1),  cpt );
+                                         }else{
+                                                 fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, k,kl+1+1,nbcode[Tvaraff[k]][lv]);
+                                                 kl++;
+                                         }
+                                 } /* end covariate */
+                         }
+                         fprintf(ficgp,"\nset out \n");
      } /* k1 */
    } /* cpt */
    /*2 eme*/
    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 */
+                                 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,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(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 */
+                                 /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-         vlv= nbcode[Tvaraff[lv]][lv];
+                                 vlv= nbcode[Tvaraff[k]][lv];
-         fprintf(ficgp," V%d=%d ",k,vlv);
+                                 fprintf(ficgp," V%d=%d ",k,vlv);
        }
        fprintf(ficgp,"\n#\n");
-     fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
+                         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*/
+                         for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
-       if(vpopbased==0)
+                                 if(vpopbased==0)
-         fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
+                                         fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
-       else
+                                 else
-         fprintf(ficgp,"\nreplot ");
+                                         fprintf(ficgp,"\nreplot ");
-       for (i=1; i<= nlstate+1 ; i ++) {
+                                 for (i=1; i<= nlstate+1 ; i ++) {
-         k=2*i;
+                                         k=2*i;
-         fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
+                                         fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
-         for (j=1; j<= nlstate+1 ; j ++) {
+                                         for (j=1; j<= nlstate+1 ; j ++) {
-           if (j==i) fprintf(ficgp," %%lf (%%lf)");
+                                                 if (j==i) fprintf(ficgp," %%lf (%%lf)");
-           else fprintf(ficgp," %%*lf (%%*lf)");
+                                                 else fprintf(ficgp," %%*lf (%%*lf)");
-         }
+                                         }
-         if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
+                                         if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
-         else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
+                                         else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
-         fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+                                         fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
-         for (j=1; j<= nlstate+1 ; j ++) {
+                                         for (j=1; j<= nlstate+1 ; j ++) {
-           if (j==i) fprintf(ficgp," %%lf (%%lf)");
+                                                 if (j==i) fprintf(ficgp," %%lf (%%lf)");
-           else fprintf(ficgp," %%*lf (%%*lf)");
+                                                 else fprintf(ficgp," %%*lf (%%*lf)");
-         }
+                                         }
-         fprintf(ficgp,"\" t\"\" w l lt 0,");
+                                         fprintf(ficgp,"\" t\"\" w l lt 0,");
-         fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+                                         fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
-         for (j=1; j<= nlstate+1 ; j ++) {
+                                         for (j=1; j<= nlstate+1 ; j ++) {
-           if (j==i) fprintf(ficgp," %%lf (%%lf)");
+                                                 if (j==i) fprintf(ficgp," %%lf (%%lf)");
-           else fprintf(ficgp," %%*lf (%%*lf)");
+                                                 else fprintf(ficgp," %%*lf (%%*lf)");
-         }
+                                         }
-         if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
+                                         if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
-         else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
+                                         else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
-       } /* state */
+                                 } /* state */
-     } /* vpopbased */
+                         } /* vpopbased */
-     fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
+                         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 */
+                                 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,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(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 */
+                                 /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
-         vlv= nbcode[Tvaraff[lv]][lv];
+                                 vlv= nbcode[Tvaraff[k]][lv];
-         fprintf(ficgp," V%d=%d ",k,vlv);
+                                 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);
        fprintf(ficgp,"set ter svg size 640, 480\n\
  plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt);
        /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
-         for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+                                 for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
-         fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+                                 fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
-         fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
+                                 fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
-         for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+                                 for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
-         fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+                                 fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
        */
        for (i=1; i< nlstate ; i ++) {
-         fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
+                                 fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
-         /*      fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
+                                 /*      fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
        }
        fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
      }
- Line 5319  plot [%.f:%.f] \"%s\" every :::%d::%d u
+ Line 5835  plot [%.f:%.f] \"%s\" every :::%d::%d u
          /* 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];
+         vlv= nbcode[Tvaraff[k]][lv];
          fprintf(ficgp," V%d=%d ",k,vlv);
        }
        fprintf(ficgp,"\n#\n");
- Line 5331  unset log y\n\
+ Line 5847  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 5354  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 5871  plot [%.f:%.f]  ", ageminpar, agemaxpar)
          /* 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];
+         vlv= nbcode[Tvaraff[k]][lv];
          fprintf(ficgp," V%d=%d ",k,vlv);
        }
        fprintf(ficgp,"\n#\n");
- Line 5398  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 5915  plot [%.f:%.f]  ", ageminpar, agemaxpar)
          /* 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];
+         vlv= nbcode[Tvaraff[k]][lv];
          fprintf(ficgp," V%d=%d ",k,vlv);
        }
        fprintf(ficgp,"\n#\n");
- Line 5423  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 5940  plot [%.f:%.f]  ", ageminpar, agemaxpar)
        fprintf(ficgp,"\nset out\n");
      } /* end cpt state*/
    } /* end covariate */
+   if(backcast == 1){
-   if(prevfcast==1){
+     /* CV back preval stable (period) for each covariate */
-   /* 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);
+         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 correspondig covariate value  */
+         for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
-           lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+           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];
+           vlv= nbcode[Tvaraff[k]][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,"PB_"),cpt,k1);
-         fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
+         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
-         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
+ set ter svg size 640, 480\n                                             \
- set ter svg size 640, 480\n\
+ unset log y\n                                                           \
- unset log y\n\
  plot [%.f:%.f]  ", ageminpar, agemaxpar);
-         for (i=1; i<= nlstate+1 ; i ++){  /* nlstate +1 p11 p21 p.1 */
+         k=3; /* Offset */
-           /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+         for (i=1; i<= nlstate ; i ++){
-           /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */
+           if(i==1)
-           /*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+             fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));
-           /*#   1       2   3    4    5      6  7   8   9   10   11 12  13   14  15 */
+           else
-           if(i==1){
+             fprintf(ficgp,", '' ");
-             fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_"));
+           /* l=(nlstate+ndeath)*(i-1)+1; */
-           }else{
+           l=(nlstate+ndeath)*(cpt-1)+1;
-             fprintf(ficgp,",\\\n '' ");
+           /* 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 *\/ */
-           if(cptcoveff ==0){ /* No covariate */
+           fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+(cpt-1)+i-1); /* a vérifier */
-             fprintf(ficgp," u 2:("); /* Age is in 2 */
+           /* for (j=2; j<= nlstate ; j ++) */
-             /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
+           /*    fprintf(ficgp,"+$%d",k+l+j-1); */
-             /*#   1       2   3   4   5  6    7  8   9   10  11  12  13  14  15  16  17  18 */
+           /*    /\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */
-             if(i==nlstate+1)
+           fprintf(ficgp,") t \"bprev(%d,%d)\" w l",i,cpt);
-               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 */
+   } /* End if backcast */
+   if(prevfcast==1){
-   /* proba elementaires */
+     /* Projection from cross-sectional to stable (period) for each covariate */
-   fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
+     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[k]][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 */
+                                                 ioffset=2; /* Age is in 2 */
+                                                 /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
+                                                 /*#   1       2   3   4   5  6    7  8   9   10  11  12  13  14  15  16  17  18 */
+                                                 /*# V1  = 1 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 */
+                                                 fprintf(ficgp," u %d:(", ioffset);
+                                                 if(i==nlstate+1)
+                                                         fprintf(ficgp," $%d/(1.-$%d)) t 'pw.%d' with line ",                    \
+                                                                                         ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+                                                 else
+                                                         fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ",                    \
+                                                                                         ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+                                         }else{ /* more than 2 covariates */
+                                                 if(cptcoveff ==1){
+                                                         ioffset=4; /* Age is in 4 */
+                                                 }else{
+                                                         ioffset=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 */
+                                                 }
+                                                 fprintf(ficgp," u %d:((",ioffset);
+                                                 kl=0;
+                                                 for (k=1; k<=cptcoveff; k++){    /* For each covariate  */
+                                                         lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to combination k1 and covariate k */
+                                                         /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
+                                                         /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
+                                                         /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
+                                                         vlv= nbcode[Tvaraff[k]][lv];
+                                                         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){
+                                                                         if(cptcoveff ==1){
+                                                                         fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[k]][lv], \
+                                                                                                         ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+                                                                         }else{
+                                                                         fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[k]][lv], \
+                                                                                                         ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+                                                                         }
+                                                                 }else{
+                                                                         if(cptcoveff ==1){
+                                                                                 fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[k]][lv], \
+                                                                                                                 ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,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[k]][lv], \
+                                                                                                                 ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+                                                                         }
+                                                                 }
+                                                         }else{ /* k < cptcoveff */
+                                                                 fprintf(ficgp,"$%d==%d && $%d==%d && ",kl, k,kl+1,nbcode[Tvaraff[k]][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");
    for(i=1,jk=1; i <=nlstate; i++){
      fprintf(ficgp,"# initial state %d\n",i);
      for(k=1; k <=(nlstate+ndeath); k++){
- Line 5591  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6167  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 5636  plot [%.f:%.f]  ", ageminpar, agemaxpar)
+ Line 6221  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 sum=0.;
+   double age;
+   double *sumnewp, *sumnewm;
+   double *agemingood, *agemaxgood; /* Currently identical for all covariates */
    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 */
+   sumnewp = vector(1,modcovmax);
+   sumnewm = vector(1,modcovmax);
+   agemingood = vector(1,modcovmax);
+   agemaxgood = vector(1,modcovmax);
+   for (cptcod=1;cptcod<=modcovmax;cptcod++){
+                 sumnewm[cptcod]=0.;
+                 sumnewp[cptcod]=0.;
+                 agemingood[cptcod]=0;
+                 agemaxgood[cptcod]=0;
+         }
+   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;
      for (age=bage; age<=fage; age++)
        for (i=1; i<=nlstate;i++)
-         for (cptcod=1;cptcod<=modcovmax;cptcod++)
+                                 for (cptcod=1;cptcod<=modcovmax;cptcod++)
-           mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
+                                         mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
      /* We keep the original values on the extreme ages bage, fage and for
         fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2
         we use a 5 terms etc. until the borders are no more concerned.
      */
      for (mob=3;mob <=mobilavrange;mob=mob+2){
        for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){
-         for (i=1; i<=nlstate;i++){
+                                 for (i=1; i<=nlstate;i++){
-           for (cptcod=1;cptcod<=modcovmax;cptcod++){
+                                         for (cptcod=1;cptcod<=modcovmax;cptcod++){
-             mobaverage[(int)age][i][cptcod] =probs[(int)age][i][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;
+                                                 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 */
+     sum=0.;
+     for (i=1; i<=nlstate;i++){
+       sum+=mobaverage[(int)agemingood[cptcod]][i][cptcod];
+     }
+     if(fabs(sum - 1.) > 1.e-3) { /* bad */
+       printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any descending age!\n",cptcod);
+       /* for (i=1; i<=nlstate;i++){ */
+       /*   mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */
+       /* } /\* i *\/ */
+     } /* end bad */
+     /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */
+                 /* 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 */
+                 sum=0.;
+                 for (i=1; i<=nlstate;i++){
+                         sum+=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
+                 }
+                 if(fabs(sum - 1.) > 1.e-3) { /* bad */
+                         printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any ascending age!\n",cptcod);
+                         /* for (i=1; i<=nlstate;i++){ */
+                         /*   mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */
+                         /* } /\* i *\/ */
+                 } /* end bad */
+                 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 5689  void prevforecast(char fileres[], double
+ Line 6382  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;
- Line 5699  void prevforecast(char fileres[], double
+ Line 6392  void prevforecast(char fileres[], double
    */
    /* 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(fileresf,"F_");
    strcat(fileresf,fileresu);
- Line 5712  void prevforecast(char fileres[], double
+ Line 6404  void prevforecast(char fileres[], double
    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 5751  void prevforecast(char fileres[], double
+ Line 6436  void prevforecast(char fileres[], double
        k=k+1;
        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",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+                                 fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
        }
        fprintf(ficresf," yearproj age");
        for(j=1; j<=nlstate+ndeath;j++){
-         for(i=1; i<=nlstate;i++)
+                                 for(i=1; i<=nlstate;i++)
            fprintf(ficresf," p%d%d",i,j);
-         fprintf(ficresf," p.%d",j);
+                                 fprintf(ficresf," wp.%d",j);
        }
        for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
-         fprintf(ficresf,"\n");
+                                 fprintf(ficresf,"\n");
-         fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);
+                                 fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);
+                                 for (agec=fage; agec>=(ageminpar-1); agec--){
-         for (agec=fage; agec>=(ageminpar-1); agec--){
+                                         nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
-           nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
+                                         nhstepm = nhstepm/hstepm;
-           nhstepm = nhstepm/hstepm;
+                                         p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
-           p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+                                         oldm=oldms;savm=savms;
-           oldm=oldms;savm=savms;
+                                         hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
-           hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
+                                         for (h=0; h<=nhstepm; h++){
-           for (h=0; h<=nhstepm; h++){
+                                                 if (h*hstepm/YEARM*stepm ==yearp) {
-             if (h*hstepm/YEARM*stepm ==yearp) {
                fprintf(ficresf,"\n");
                for(j=1;j<=cptcoveff;j++)
                  fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-               fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
+                                                         fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
-             }
+                                                 }
-             for(j=1; j<=nlstate+ndeath;j++) {
+                                                 for(j=1; j<=nlstate+ndeath;j++) {
-               ppij=0.;
+                                                         ppij=0.;
-               for(i=1; i<=nlstate;i++) {
+                                                         for(i=1; i<=nlstate;i++) {
-                 if (mobilav==1)
+                                                                 if (mobilav==1)
-                   ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod];
+                                                                         ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod];
-                 else {
+                                                                 else {
-                   ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];
+                                                                         ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];
-                 }
+                                                                 }
-                 if (h*hstepm/YEARM*stepm== yearp) {
+                                                                 if (h*hstepm/YEARM*stepm== yearp) {
-                   fprintf(ficresf," %.3f", p3mat[i][j][h]);
+                                                                         fprintf(ficresf," %.3f", p3mat[i][j][h]);
-                 }
+                                                                 }
-               } /* end i */
+                                                         } /* end i */
-               if (h*hstepm/YEARM*stepm==yearp) {
+                                                         if (h*hstepm/YEARM*stepm==yearp) {
-                 fprintf(ficresf," %.3f", ppij);
+                                                                 fprintf(ficresf," %.3f", ppij);
-               }
+                                                         }
-             }/* end j */
+                                                 }/* end j */
-           } /* end h */
+                                         } /* end h */
-           free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+                                         free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
-         } /* end agec */
+                                 } /* end agec */
        } /* end yearp */
      } /* 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 5816  void populforecast(char fileres[], doubl
+ Line 6628  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 5838  void populforecast(char fileres[], doubl
+ Line 6650  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 5853  void populforecast(char fileres[], doubl
+ Line 6665  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 5865  void populforecast(char fileres[], doubl
+ Line 6677  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 5885  void populforecast(char fileres[], doubl
+ Line 6697  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 5908  void populforecast(char fileres[], doubl
+ Line 6720  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 5953  void populforecast(char fileres[], doubl
+ Line 6766  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 5967  void populforecast(char fileres[], doubl
+ Line 6780  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 6242  int readdata(char datafile[], int firsto
+ Line 7055  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 6620  int decodemodel ( char model[], int last
+ Line 7433  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 6916  void syscompilerinfo(int logged)
+ Line 7734  void syscompilerinfo(int logged)
  #endif
-  }
+ }
-  int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
+ int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
    /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
    int i, j, k, i1 ;
    /* double ftolpl = 1.e-10; */
- Line 6939  void syscompilerinfo(int logged)
+ Line 7757  void syscompilerinfo(int logged)
    for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
    fprintf(ficrespl,"\n");
-     /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
+   /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
-     agebase=ageminpar;
+   agebase=ageminpar;
-     agelim=agemaxpar;
+   agelim=agemaxpar;
-     i1=pow(2,cptcoveff);
+   i1=pow(2,cptcoveff);
-     if (cptcovn < 1){i1=1;}
+   if (cptcovn < 1){i1=1;}
-     for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+   for(cptcov=1,k=0;cptcov<=i1;cptcov++){
      /* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
-       //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+     //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
-         k=k+1;
+     k=k+1;
-         /* to clean */
+     /* to clean */
-         //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+     //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
-         fprintf(ficrespl,"#******");
+     fprintf(ficrespl,"#******");
-         printf("#******");
+     printf("#******");
-         fprintf(ficlog,"#******");
+     fprintf(ficlog,"#******");
-         for(j=1;j<=cptcoveff;j++) {
+     for(j=1;j<=cptcoveff;j++) {
-           fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+       fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-           printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+       printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-           fprintf(ficlog," 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(ficrespl,"******\n");
+     fprintf(ficrespl,"******\n");
-         printf("******\n");
+     printf("******\n");
-         fprintf(ficlog,"******\n");
+     fprintf(ficlog,"******\n");
-         fprintf(ficrespl,"#Age ");
+     fprintf(ficrespl,"#Age ");
-         for(j=1;j<=cptcoveff;j++) {
+     for(j=1;j<=cptcoveff;j++) {
-           fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+       fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-         }
+     }
-         for(i=1; i<=nlstate;i++) fprintf(ficrespl,"  %d-%d   ",i,i);
+     for(i=1; i<=nlstate;i++) fprintf(ficrespl,"  %d-%d   ",i,i);
-         fprintf(ficrespl,"Total Years_to_converge\n");
+     fprintf(ficrespl,"Total Years_to_converge\n");
-         for (age=agebase; age<=agelim; age++){
+     for (age=agebase; age<=agelim; age++){
-         /* for (age=agebase; age<=agebase; age++){ */
+       /* for (age=agebase; age<=agebase; age++){ */
-           prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
+       prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
-           fprintf(ficrespl,"%.0f ",age );
+       fprintf(ficrespl,"%.0f ",age );
-           for(j=1;j<=cptcoveff;j++)
+       for(j=1;j<=cptcoveff;j++)
-             fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+         fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
-           tot=0.;
+       tot=0.;
-           for(i=1; i<=nlstate;i++){
+       for(i=1; i<=nlstate;i++){
-             tot +=  prlim[i][i];
+         tot +=  prlim[i][i];
-             fprintf(ficrespl," %.5f", prlim[i][i]);
+         fprintf(ficrespl," %.5f", prlim[i][i]);
-           }
+       }
-           fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
+       fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
-         } /* Age */
+     } /* Age */
-         /* was end of cptcod */
+     /* was end of cptcod */
-     } /* cptcov */
+   } /* cptcov */
-         return 0;
+   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 7016  int hPijx(double *p, int bage, int fage)
+ Line 7930  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 7057  int hPijx(double *p, int bage, int fage)
+ Line 7971  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 7091  int main(int argc, char *argv[])
+ Line 8082  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 7111  int main(int argc, char *argv[])
+ Line 8103  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 7121  int main(int argc, char *argv[])
+ Line 8114  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 7132  int main(int argc, char *argv[])
+ Line 8126  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 7568  This is probably because your covariance
+ Line 8564  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 7609  Please run with mle=-1 to get a correct
+ Line 8605  Please run with mle=-1 to get a correct
      }
      fprintf(ficres,"#%s\n",version);
    }    /* End of mle != -3 */
    /*  Main data
     */
    n= lastobs;
- Line 7784  Please run with mle=-1 to get a correct
+ Line 8780  Please run with mle=-1 to get a correct
       *              bbbbbbbb
       *              76543210
       *   h-1        00000101 (6-1=5)
-      *(h-1)>>(k-1)= 00000001 >> (2-1) = 1 right shift
+      *(h-1)>>(k-1)= 00000010 >> (2-1) = 1 right shift
       *           &
       *     1        00000001 (1)
-      *              00000001        = 1 & ((h-1) >> (k-1))
+      *              00000000        = 1 & ((h-1) >> (k-1))
-      *          +1= 00000010 =2
+      *          +1= 00000001 =1
       *
       * h=14, k=3 => h'=h-1=13, k'=k-1=2
       *          h'      1101 =2^3+2^2+0x2^1+2^0
- Line 7933  Interval (in months) between two waves:
+ Line 8929  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 7952  Interval (in months) between two waves:
+ Line 8947  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 8367  Please run with mle=-1 to get a correct
+ Line 9362  Please run with mle=-1 to get a correct
      if((num_filled=sscanf(line,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm, &ftolpl)) !=EOF){
      if (num_filled != 6) {
-       printf("Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n");
+       printf("Error: Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n, your line=%s . Probably you are running an older format.\n",line);
-       printf("but line=%s\n",line);
+       fprintf(ficlog,"Error: Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n, your line=%s . Probably you are running an older format.\n",line);
        goto end;
      }
      printf("agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",ageminpar,agemaxpar, bage, fage, estepm, ftolpl);
- Line 8434  Please run with mle=-1 to get a correct
+ Line 9429  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);
+     fprintf(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);
+     fprintf(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);
+     fprintf(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 8448  Please run with mle=-1 to get a correct
+ Line 9456  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, prevfcast, 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,prevfcast,estepm, \
+                  model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,backcast, estepm, \
                   jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2);
     /*------------ free_vector  -------------*/
- Line 8478  Please run with mle=-1 to get a correct
+ Line 9486  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);
-   /*-------------- Variance of one-step probabilities---*/
+     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+ndeath;j++) /* ndeath is useless but a necessity to be compared with mobaverages */
-         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 ) {
+       mobaverages= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
+                         for(i=1;i<=AGESUP;i++)
+                                 for(j=1;j<=nlstate;j++)
+                                         for(k=1;k<=ncovcombmax;k++)
+                                                 mobaverages[i][j][k]=0.;
+       mobaverage=mobaverages;
+       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) --------------*/
+       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);
+     }
      /* ------ 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);
-     /*  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);
      free_imatrix(dh,1,lastpass-firstpass+2,1,imx);
      free_imatrix(bh,1,lastpass-firstpass+2,1,imx);
      free_imatrix(mw,1,lastpass-firstpass+2,1,imx);
-     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 8545  Please run with mle=-1 to get a correct
+ Line 9573  Please run with mle=-1 to get a correct
      }
      printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout);
      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);
-         oldm=oldms;savm=savms;
-         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);
+       eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
-       /*}*/
+       oldm=oldms;savm=savms;
+       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);
      }
      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 8577  Please run with mle=-1 to get a correct
+ Line 9602  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 8612  Please run with mle=-1 to get a correct
+ Line 9637  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 8737  Please run with mle=-1 to get a correct
+ Line 9762  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(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
+       free_ma3x(mobaverages,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,nlstate+ndeath, 1,ncovcombmax);
    }  /* mle==-3 arrives here for freeing */
   /* endfree:*/
      free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */
- Line 8784  Please run with mle=-1 to get a correct
+ Line 9810  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>