--- imach/src/imach.c	2018/04/27 14:27:04	1.286
+++ imach/src/imach.c	2019/05/09 15:17:34	1.293
@@ -1,6 +1,29 @@
-/* $Id: imach.c,v 1.286 2018/04/27 14:27:04 brouard Exp $
+/* $Id: imach.c,v 1.293 2019/05/09 15:17:34 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.293  2019/05/09 15:17:34  brouard
+  *** empty log message ***
+
+  Revision 1.292  2019/05/09 14:17:20  brouard
+  Summary: Some updates
+
+  Revision 1.291  2019/05/09 13:44:18  brouard
+  Summary: Before ncovmax
+
+  Revision 1.290  2019/05/09 13:39:37  brouard
+  Summary: 0.99r18 unlimited number of individuals
+
+  The number n which was limited to 20,000 cases is now unlimited, from firstobs to lastobs. If the number is too for the virtual memory, probably an error will occur.
+
+  Revision 1.289  2018/12/13 09:16:26  brouard
+  Summary: Bug for young ages (<-30) will be in r17
+
+  Revision 1.288  2018/05/02 20:58:27  brouard
+  Summary: Some bugs fixed
+
+  Revision 1.287  2018/05/01 17:57:25  brouard
+  Summary: Bug fixed by providing frequencies only for non missing covariates
+
   Revision 1.286  2018/04/27 14:27:04  brouard
   Summary: some minor bugs
 
@@ -1040,14 +1063,15 @@ typedef struct {
 #define NINTERVMAX 8
 #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */
 #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
-#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */
+#define NCOVMAX 20  /**< Maximum number of covariates, including generated covariates V1*V2 */
 #define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
 /*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/
 #define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 
-#define MAXN 20000
+/*#define MAXN 20000 */ /* Should by replaced by nobs, real number of observations and unlimited */
 #define YEARM 12. /**< Number of months per year */
 /* #define AGESUP 130 */
-#define AGESUP 150
+/* #define AGESUP 150 */
+#define AGESUP 200
 #define AGEINF 0
 #define AGEMARGE 25 /* Marge for agemin and agemax for(iage=agemin-AGEMARGE; iage <= agemax+3+AGEMARGE; iage++) */
 #define AGEBASE 40
@@ -1063,12 +1087,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.286 2018/04/27 14:27:04 brouard Exp $ */
+/* $Id: imach.c,v 1.293 2019/05/09 15:17:34 brouard Exp $ */
 /* $State: Exp $ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
 char copyright[]="April 2018,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018";
-char fullversion[]="$Revision: 1.286 $ $Date: 2018/04/27 14:27:04 $"; 
+char fullversion[]="$Revision: 1.293 $ $Date: 2019/05/09 15:17:34 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1091,6 +1115,7 @@ int nqfveff=0; /**< nqfveff Number of Qu
 int ntveff=0; /**< ntveff number of effective time varying variables */
 int nqtveff=0; /**< ntqveff number of effective time varying quantitative variables */
 int cptcov=0; /* Working variable */
+int nobs=10;  /* Number of observations in the data lastobs-firstobs */
 int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */
 int npar=NPARMAX;
 int nlstate=2; /* Number of live states */
@@ -2574,6 +2599,7 @@ void powell(double p[], double **xi, int
   double **newm;
   double agefin, delaymax=200. ; /* 100 Max number of years to converge */
   int ncvloop=0;
+  int first=0;
   
   min=vector(1,nlstate);
   max=vector(1,nlstate);
@@ -2670,10 +2696,14 @@ void powell(double p[], double **xi, int
       free_vector(meandiff,1,nlstate);
       return prlim;
     }
-  } /* age loop */
+  } /* agefin loop */
     /* After some age loop it doesn't converge */
-  printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
-Earliest 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);
+  if(!first){
+    first=1;
+    printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d). Others in log file only...\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM),  (int)(age-stepm/YEARM), (int)delaymax);
+  }
+  fprintf(ficlog, "Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d).\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM),  (int)(age-stepm/YEARM), (int)delaymax);
+
   /* 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);
@@ -2739,7 +2769,8 @@ Earliest age to start was %d-%d=%d, ncvl
   /* 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 */
+  /* for(agefin=age; agefin<AGESUP; agefin=agefin+stepm/YEARM){ /\* A changer en age *\/ */
+  for(agefin=age; agefin<FMIN(AGESUP,age+delaymax); 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 */
@@ -2853,9 +2884,9 @@ Earliest age to start was %d-%d=%d, ncvl
       free_vector(meandiff,1,nlstate);
       return bprlim;
     }
-  } /* age loop */
+  } /* agefin loop */
     /* After some age loop it doesn't converge */
-  if(first){
+  if(!first){
     first=1;
     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'. Others in log file only...\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);
@@ -2938,8 +2969,24 @@ double **pmij(double **ps, double *cov,
       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]); */
@@ -2967,7 +3014,7 @@ double **pmij(double **ps, double *cov,
   double **out, **pmij();
   double sumnew=0.;
   double agefin;
-  double k3=0.; /* constant of the w_x diagonal matrixe (in order for B to sum to 1 even for death state) */
+  double k3=0.; /* constant of the w_x diagonal matrix (in order for B to sum to 1 even for death state) */
   double **dnewm, **dsavm, **doldm;
   double **bbmij;
   
@@ -2986,7 +3033,7 @@ double **pmij(double **ps, double *cov,
   /* outputs pmmij which is a stochastic matrix in row */
 
   /* Diag(w_x) */
-  /* Problem with prevacurrent which can be zero */
+  /* Rescaling the cross-sectional prevalence: Problem with prevacurrent which can be zero */
   sumnew=0.;
   /*for (ii=1;ii<=nlstate+ndeath;ii++){*/
   for (ii=1;ii<=nlstate;ii++){ /* Only on live states */
@@ -3013,10 +3060,10 @@ double **pmij(double **ps, double *cov,
   }
   /* End doldm, At the end doldm is diag[(w_i)] */
   
-  /* left Product of this diag matrix by pmmij=Px (dnewm=dsavm*doldm) */
-  bbmij=matprod2(dnewm, doldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, pmmij); /* Bug Valgrind */
+  /* Left product of this diag matrix by pmmij=Px (dnewm=dsavm*doldm): diag[(w_i)*Px */
+  bbmij=matprod2(dnewm, doldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, pmmij); /* was a Bug Valgrind */
 
-  /* Diag(Sum_i w^i_x p^ij_x */
+  /* Diag(Sum_i w^i_x p^ij_x, should be the prevalence at age x+stepm */
   /* w1 p11 + w2 p21 only on live states N1./N..*N11/N1. + N2./N..*N21/N2.=(N11+N21)/N..=N.1/N.. */
   for (j=1;j<=nlstate+ndeath;j++){
     sumnew=0.;
@@ -3034,7 +3081,7 @@ double **pmij(double **ps, double *cov,
     } /*End ii */
   } /* End j, At the end dsavm is diag[1/(w_1p1i+w_2 p2i)] for ALL states even if the sum is only for live states */
 
-  ps=matprod2(ps, dnewm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dsavm); /* Bug Valgrind */
+  ps=matprod2(ps, dnewm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dsavm); /* was a Bug Valgrind */
   /* ps is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */
   /* end bmij */
   return ps; /*pointer is unchanged */
@@ -3866,7 +3913,7 @@ return -l;
 
 
 /*************** function likelione ***********/
-void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*funcone)(double []))
+void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*func)(double []))
 {
   /* This routine should help understanding what is done with 
      the selection of individuals/waves and
@@ -3890,7 +3937,7 @@ void likelione(FILE *ficres,double p[],
     fprintf(ficresilk," -2*gipw/gsw*weight*ll(total)\n");
   }
 
-  *fretone=(*funcone)(p);
+  *fretone=(*func)(p);
   if(*globpri !=0){
     fclose(ficresilk);
     if (mle ==0)
@@ -5052,7 +5099,7 @@ void prevalence(double ***probs, double
   /*j=cptcoveff;*/
   if (cptcovn<1) {j=1;ncodemax[1]=1;}
   
-  first=1;
+  first=0;
   for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */
     for (i=1; i<=nlstate; i++)  
       for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++)
@@ -5110,12 +5157,11 @@ void prevalence(double ***probs, double
 	  if(posprop>=1.e-5){ 
 	    probs[i][jk][j1]= prop[jk][i]/posprop;
 	  } else{
-	    if(first==1){
-	      first=0;
+	    if(!first){
+	      first=1;
 	      printf("Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]);
-	      fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]);
 	    }else{
-	      fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]);
+	      fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases.\n",jk,i,jk, j1,probs[i][jk][j1]);
 	    }
 	  }
 	} 
@@ -5366,7 +5412,6 @@ void  concatwav(int wav[], int **dh, int
        nbcode[k][j]=0; /* Valgrind */
 
    /* Loop on covariates without age and products and no quantitative variable */
-   /* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */
    for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */
      for (j=-1; (j < maxncov); j++) Ndum[j]=0;
      if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ 
@@ -5384,7 +5429,11 @@ void  concatwav(int wav[], int **dh, int
 	     modmaxcovj=ij; 
 	   else if (ij < modmincovj) 
 	     modmincovj=ij; 
-	   if ((ij < -1) && (ij > NCOVMAX)){
+	   if (ij <0 || ij >1 ){
+	     printf("Information, IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i);
+	     fprintf(ficlog,"Information, currently IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i);
+	   }
+	   if ((ij < -1) || (ij > NCOVMAX)){
 	     printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
 	     exit(1);
 	   }else
@@ -5430,12 +5479,18 @@ void  concatwav(int wav[], int **dh, int
 	 /* nbcode[Tvar[j]][3]=2; */
 	 /* To be continued (not working yet). */
 	 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*/
+
+	 /* 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*/
+	 /* Skipping the case of missing values by reducing nbcode to 0 and 1 and not -1, 0, 1 */
+	 /* model=V1+V2+V3, if V2=-1, 0 or 1, then nbcode[2][1]=0 and nbcode[2][2]=1 instead of
+	  * nbcode[2][1]=-1, nbcode[2][2]=0 and nbcode[2][3]=1 */
+	 /*, could be restored in the future */
+	 for (i=0; i<=1; 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;
 	   }
 	   ij++;
-	   nbcode[Tvar[k]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/
+	   nbcode[Tvar[k]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1 . Could be -1*/
 	   cptcode = ij; /* New max modality for covar j */
 	 } /* end of loop on modality i=-1 to 1 or more */
 	 break;
@@ -5451,21 +5506,7 @@ void  concatwav(int wav[], int **dh, int
 	 break;
        } /* end switch */
      } /* end dummy test */
-    
-     /*   for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */
-     /* 	/\*recode from 0 *\/ */
-     /* 				     k is a modality. If we have model=V1+V1*sex  */
-     /* 				     then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */
-     /* 				  But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */
-     /* 	} */
-     /* 	/\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */
-     /* 	if (ij > ncodemax[j]) { */
-     /* 	  printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]);  */
-     /* 	  fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */
-     /* 	  break; */
-     /* 	} */
-     /*   }  /\* 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*/  
+   } /* end of loop on model-covariate k. nbcode[Tvark][1]=-1, nbcode[Tvark][1]=0 and nbcode[Tvark][2]=1 sets the value of covariate k*/  
   
    for (k=-1; k< maxncov; k++) Ndum[k]=0; 
    /* Look at fixed dummy (single or product) covariates to check empty modalities */
@@ -5849,6 +5890,7 @@ void  concatwav(int wav[], int **dh, int
    double **dnewm,**doldm;
    double **dnewmp,**doldmp;
    int i, j, nhstepm, hstepm, h, nstepm ;
+   int first=0;
    int k;
    double *xp;
    double **gp, **gm;  /**< for var eij */
@@ -5973,7 +6015,7 @@ void  concatwav(int wav[], int **dh, int
        }
        /**< Computes the prevalence limit with parameter theta shifted of delta up to ftolpl precision and 
 	* returns into prlim .
-	*/		
+	*/
        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres);
 
        /* If popbased = 1 we use crossection prevalences. Previous step is useless but prlim is created */
@@ -5989,7 +6031,7 @@ void  concatwav(int wav[], int **dh, int
        /**< Computes the shifted transition matrix \f$ {}{h}_p^{ij}_x\f$ at horizon h.
 	*/			
        hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres);  /* Returns p3mat[i][j][h] for h=0 to nhstepm */
-       /**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}_x\f$, which are the probability
+       /**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}x\f$, which are the probability
 	* at horizon h in state j including mortality.
 	*/
        for(j=1; j<= nlstate; j++){
@@ -6011,7 +6053,7 @@ void  concatwav(int wav[], int **dh, int
 			
        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, nres);
 			
        if (popbased==1) {
@@ -6197,7 +6239,7 @@ void  concatwav(int wav[], int **dh, int
   int theta;
   
   pstamp(ficresvpl);
-  fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n");
+  fprintf(ficresvpl,"# Standard deviation of period (forward stable) prevalences \n");
   fprintf(ficresvpl,"# Age ");
   if(nresult >=1)
     fprintf(ficresvpl," Result# ");
@@ -6226,20 +6268,20 @@ void  concatwav(int wav[], int **dh, int
       for(i=1; i<=npar; i++){ /* Computes gradient */
 	xp[i] = x[i] + (i==theta ?delti[theta]:0);
       }
-      if((int)age==79 ||(int)age== 80 ||(int)age== 81 )
-	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres);
-      else
-	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres);
+      /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) */
+      /* 	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); */
+      /* else */
+      prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres);
       for(i=1;i<=nlstate;i++){
 	gp[i] = prlim[i][i];
 	mgp[theta][i] = prlim[i][i];
       }
       for(i=1; i<=npar; i++) /* Computes gradient */
 	xp[i] = x[i] - (i==theta ?delti[theta]:0);
-      if((int)age==79 ||(int)age== 80 ||(int)age== 81 )
-	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres);
-      else
-	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres);
+      /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) */
+      /* 	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); */
+      /* else */
+      prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres);
       for(i=1;i<=nlstate;i++){
 	gm[i] = prlim[i][i];
 	mgm[theta][i] = prlim[i][i];
@@ -6288,8 +6330,11 @@ void  concatwav(int wav[], int **dh, int
     fprintf(ficresvpl,"%.0f ",age );
     if(nresult >=1)
       fprintf(ficresvpl,"%d ",nres );
-    for(i=1; i<=nlstate;i++)
+    for(i=1; i<=nlstate;i++){
       fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age]));
+      /* for(j=1;j<=nlstate;j++) */
+      /* 	fprintf(ficresvpl," %d %.5f ",j,prlim[j][i]); */
+    }
     fprintf(ficresvpl,"\n");
     free_vector(gp,1,nlstate);
     free_vector(gm,1,nlstate);
@@ -6506,7 +6551,7 @@ void varprob(char optionfilefiname[], do
    fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
    fprintf(fichtm,"\n");
 
-   fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back. %s</li>\n",optionfilehtmcov,optionfilehtmcov);
+   fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back. File %s</li>\n",optionfilehtmcov,optionfilehtmcov);
    fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov);
    fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \
 and drawn. It helps understanding how is the covariance between two incidences.\
@@ -6799,10 +6844,10 @@ void printinghtml(char fileresu[], char
  - 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",
+ - Period (forward) 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",
+ - Backward 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): \
@@ -6908,31 +6953,31 @@ divided by h: <sub>h</sub>P<sub>ij</sub>
 <img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); 
      /* Survival functions (period) in state j */
      for(cpt=1; cpt<=nlstate;cpt++){
-       fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+       fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
 <img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
      }
      /* State specific survival functions (period) */
      for(cpt=1; cpt<=nlstate;cpt++){
-       fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\
- Or probability to survive in various states (1 to %d) being in state %d at different ages.	\
+       fprintf(fichtm,"<br>\n- Survival functions in state %d and in any other live state (total).\
+ And probability to be observed in various states (up to %d) being in state %d at different ages.	\
  <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> <img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
      }
-     /* Period (stable) prevalence in each health state */
+     /* Period (forward 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 for a person being in state (1 to %d) at different ages, to be in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
 <img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
      }
      if(backcast==1){
-       /* Period (stable) back prevalence in each health state */
+       /* Backward prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){
 	 fprintf(fichtm,"<br>\n- Convergence to mixed (stable) back prevalence in state %d. Or probability for a person to be in state %d at a younger age, knowing that she/he was in state (1 to %d) at different older ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
 <img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
        }
      }
      if(prevfcast==1){
-       /* Projection of prevalence up to period (stable) prevalence in each health state */
+       /* Projection of prevalence up to period (forward stable) prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){
-	 fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+	 fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) forward prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
 <img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, dateproj1, dateproj2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
        }
      }
@@ -6985,13 +7030,13 @@ See page 'Matrix of variance-covariance
    <a href=\"%s\">%s</a> <br>\n</li>",
 	   estepm,subdirf2(fileresu,"STDE_"),subdirf2(fileresu,"STDE_"));
    fprintf(fichtm,"\
- - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), e<sup>ij</sup> are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a><br>\n",
+ - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), e<sup>ij</sup> are weighted by the forward (period) prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a><br>\n",
 	   estepm, subdirf2(fileresu,"V_"),subdirf2(fileresu,"V_"));
    fprintf(fichtm,"\
  - Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a> <br>\n",
 	   estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"T_"));
    fprintf(fichtm,"\
- - Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\
+ - Standard deviation of forward (period) prevalences: <a href=\"%s\">%s</a> <br>\n",\
 	   subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_"));
 
 /*  if(popforecast==1) fprintf(fichtm,"\n */
@@ -7127,7 +7172,7 @@ void printinggnuplot(char fileresu[], ch
 	  continue;
 	/* We are interested in selected combination by the resultline */
 	/* printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); */
-	fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files  and live state =%d ", cpt);
+	fprintf(ficgp,"\n# 1st: Forward (stable period) prevalence with CI: 'VPL_' files  and live state =%d ", cpt);
 	strcpy(gplotlabel,"(");
 	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 value of the covariate corresponding to k1 combination */
@@ -7165,7 +7210,7 @@ void printinggnuplot(char fileresu[], ch
 	  if (i==cpt) 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==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres);
+	fprintf(ficgp,"\" t\"Forward prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres);
 	for (i=1; i<= nlstate ; i ++) {
 	  if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
 	  else fprintf(ficgp," %%*lf (%%*lf)");
@@ -7494,7 +7539,7 @@ set ter svg size 640, 480\nunset log y\n
       continue;
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state of arrival */
       strcpy(gplotlabel,"(");      
-      fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
+      fprintf(ficgp,"\n#\n#\n#CV preval stable (forward): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
 	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
 	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
@@ -7538,14 +7583,14 @@ set ter svg size 640, 480\nunset log y\n
   
 /* 7eme */
   if(backcast == 1){
-    /* CV back preval stable (period) for each covariate */
+    /* CV backward prevalence  for each covariate */
     for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
       if(m != 1 && TKresult[nres]!= k1)
 	continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life origin state */
 	strcpy(gplotlabel,"(");      
- 	fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pijb' files, covariatecombination#=%d state=%d",k1, cpt);
+ 	fprintf(ficgp,"\n#\n#\n#CV Backward stable prevalence: 'pijb' files, covariatecombination#=%d state=%d",k1, cpt);
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
 	  lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
 	  /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
@@ -7593,7 +7638,7 @@ set ter svg size 640, 480\nunset log y\n
   
   /* 8eme */
   if(prevfcast==1){
-    /* Projection from cross-sectional to stable (period) for each covariate */
+    /* Projection from cross-sectional to forward stable (period) prevalence for each covariate */
     
     for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
@@ -7601,7 +7646,7 @@ set ter svg size 640, 480\nunset log y\n
 	continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
 	strcpy(gplotlabel,"(");      
-	fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
+	fprintf(ficgp,"\n#\n#\n#Projection of prevalence to forward stable prevalence (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 */
@@ -8035,6 +8080,7 @@ set ter svg size 640, 480\nunset log y\n
    int modcovmax =1;
    int mobilavrange, mob;
    int iage=0;
+   int firstA1=0, firstA2=0;
 
    double sum=0., sumr=0.;
    double age;
@@ -8132,6 +8178,7 @@ set ter svg size 640, 480\nunset log y\n
      } /* age */
      /* Thus we have agemingood and agemaxgood as well as goodr for raw (preobs) */
      /* but they will change */
+     firstA1=0;firstA2=0;
      for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, filling up to the youngest */
        sumnewm[cptcod]=0.;
        sumnewmr[cptcod]=0.;
@@ -8164,11 +8211,19 @@ set ter svg size 640, 480\nunset log y\n
 	 sumr+=probs[(int)age][i][cptcod];
        }
        if(fabs(sum - 1.) > 1.e-3) { /* bad */
-	 printf("Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage);
+	 if(!firstA1){
+	   firstA1=1;
+	   printf("Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d. Others in log file...\n",cptcod,sumr, (int)age, (int)bage);
+	 }
+	 fprintf(ficlog,"Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage);
        } /* end bad */
        /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */
        if(fabs(sumr - 1.) > 1.e-3) { /* bad */
-	 printf("Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage);
+	 if(!firstA2){
+	   firstA2=1;
+	   printf("Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d. Others in log file...\n",cptcod,sumr, (int)age, (int)bage);
+	 }
+	 fprintf(ficlog,"Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage);
        } /* end bad */
      }/* age */
 
@@ -8543,7 +8598,7 @@ set ter svg size 640, 480\nunset log y\n
 
 /* Variance of prevalence limit: varprlim */
  void varprlim(char fileresu[], int nresult, double ***prevacurrent, int mobilavproj, double bage, double fage, double **prlim, int *ncvyearp, double ftolpl, double p[], double **matcov, double *delti, int stepm, int cptcoveff){
-    /*------- Variance of period (stable) prevalence------*/   
+    /*------- Variance of forward period (stable) prevalence------*/   
  
    char fileresvpl[FILENAMELENGTH];  
    FILE *ficresvpl;
@@ -8554,11 +8609,11 @@ set ter svg size 640, 480\nunset log y\n
     strcpy(fileresvpl,"VPL_");
     strcat(fileresvpl,fileresu);
     if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
-      printf("Problem with variance of period (stable) prevalence  resultfile: %s\n", fileresvpl);
+      printf("Problem with variance of forward period (stable) prevalence  resultfile: %s\n", fileresvpl);
       exit(0);
     }
-    printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout);
-    fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog);
+    printf("Computing Variance-covariance of forward period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout);
+    fprintf(ficlog, "Computing Variance-covariance of forward period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog);
     
     /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
@@ -8595,8 +8650,8 @@ set ter svg size 640, 480\nunset log y\n
     }
     
     fclose(ficresvpl);
-    printf("done variance-covariance of period prevalence\n");fflush(stdout);
-    fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);
+    printf("done variance-covariance of forward period prevalence\n");fflush(stdout);
+    fprintf(ficlog,"done variance-covariance of forward period prevalence\n");fflush(ficlog);
 
  }
 /* Variance of back prevalence: varbprlim */
@@ -10166,8 +10221,10 @@ BOOL IsWow64()
 #endif
 
 void syscompilerinfo(int logged)
- {
-   /* #include "syscompilerinfo.h"*/
+{
+#include <stdint.h>
+
+  /* #include "syscompilerinfo.h"*/
    /* command line Intel compiler 32bit windows, XP compatible:*/
    /* /GS /W3 /Gy
       /Zc:wchar_t /Zi /O2 /Fd"Release\vc120.pdb" /D "WIN32" /D "NDEBUG" /D
@@ -10202,6 +10259,8 @@ void syscompilerinfo(int logged)
       /ManifestFile:"Release\IMaCh.exe.intermediate.manifest" /OPT:ICF
       /NOLOGO /TLBID:1
    */
+
+
 #if defined __INTEL_COMPILER
 #if defined(__GNUC__)
 	struct utsname sysInfo;  /* For Intel on Linux and OS/X */
@@ -10218,8 +10277,6 @@ void syscompilerinfo(int logged)
    }
 #endif
 
-#include <stdint.h>
-
    printf("Compiled with:");if(logged)fprintf(ficlog,"Compiled with:");
 #if defined(__clang__)
    printf(" Clang/LLVM");if(logged)fprintf(ficlog," Clang/LLVM");	/* Clang/LLVM. ---------------------------------------------- */
@@ -10326,7 +10383,7 @@ void syscompilerinfo(int logged)
 }
 
 int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
-  /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
+  /*--------------- Prevalence limit  (forward period or forward stable prevalence) --------------*/
   int i, j, k, i1, k4=0, nres=0 ;
   /* double ftolpl = 1.e-10; */
   double age, agebase, agelim;
@@ -10335,13 +10392,13 @@ int prevalence_limit(double *p, double *
   strcpy(filerespl,"PL_");
   strcat(filerespl,fileresu);
   if((ficrespl=fopen(filerespl,"w"))==NULL) {
-    printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
-    fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
+    printf("Problem with forward period (stable) prevalence resultfile: %s\n", filerespl);return 1;
+    fprintf(ficlog,"Problem with forward period (stable) prevalence resultfile: %s\n", filerespl);return 1;
   }
-  printf("\nComputing period (stable) prevalence: result on file '%s' \n", filerespl);
-  fprintf(ficlog,"\nComputing period (stable) prevalence: result on file '%s' \n", filerespl);
+  printf("\nComputing forward period (stable) prevalence: result on file '%s' \n", filerespl);
+  fprintf(ficlog,"\nComputing forward period (stable) prevalence: result on file '%s' \n", filerespl);
   pstamp(ficrespl);
-  fprintf(ficrespl,"# Period (stable) prevalence. Precision given by ftolpl=%g \n", ftolpl);
+  fprintf(ficrespl,"# Forward period (stable) prevalence. Precision given by ftolpl=%g \n", ftolpl);
   fprintf(ficrespl,"#Age ");
   for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
   fprintf(ficrespl,"\n");
@@ -10416,7 +10473,7 @@ int prevalence_limit(double *p, double *
 }
 
 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) --------------*/
+	/*--------------- Back Prevalence limit  (backward stable prevalence) --------------*/
 	
 	/* Computes the back prevalence limit  for any combination	of covariate values 
    * at any age between ageminpar and agemaxpar
@@ -10431,13 +10488,13 @@ int back_prevalence_limit(double *p, dou
   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("Problem with backward prevalence resultfile: %s\n", fileresplb);return 1;
+    fprintf(ficlog,"Problem with backward 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);
+  printf("Computing backward prevalence: result on file '%s' \n", fileresplb);
+  fprintf(ficlog,"Computing backward prevalence: result on file '%s' \n", fileresplb);
   pstamp(ficresplb);
-  fprintf(ficresplb,"# Period (stable) back prevalence. Precision given by ftolpl=%g \n", ftolpl);
+  fprintf(ficresplb,"# Backward 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");
@@ -10626,7 +10683,7 @@ int hPijx(double *p, int bage, int fage)
   /*if (stepm<=24) stepsize=2;*/
   
   /* agelim=AGESUP; */
-  ageminl=30;
+  ageminl=AGEINF; /* was 30 */
   hstepm=stepsize*YEARM; /* Every year of age */
   hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
   
@@ -10705,7 +10762,8 @@ int main(int argc, char *argv[])
   double ssval;
 #endif
   int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
-  int i,j, k, n=MAXN,iter=0,m,size=100, cptcod;
+  int i,j, k, iter=0,m,size=100, cptcod; /* Suppressing because nobs */
+  /* int i,j, k, n=MAXN,iter=0,m,size=100, cptcod; */
   int ncvyear=0; /* Number of years needed for the period prevalence to converge */
   int jj, ll, li, lj, lk;
   int numlinepar=0; /* Current linenumber of parameter file */
@@ -10740,7 +10798,7 @@ int main(int argc, char *argv[])
   
   char pathr[MAXLINE], pathimach[MAXLINE]; 
   char *tok, *val; /* pathtot */
-  int firstobs=1, lastobs=10;
+  int firstobs=1, lastobs=10; /* nobs = lastobs-firstobs declared globally ;*/
   int c,  h , cpt, c2;
   int jl=0;
   int i1, j1, jk, stepsize=0;
@@ -10748,7 +10806,7 @@ 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 backcast=0; /* defined as global for mlikeli and mle*/
   int mobilav=0,popforecast=0;
   int hstepm=0, nhstepm=0;
   int agemortsup;
@@ -11107,10 +11165,10 @@ int main(int argc, char *argv[])
   ungetc(c,ficpar);
 
    
-  covar=matrix(0,NCOVMAX,1,n);  /**< used in readdata */
-  if(nqv>=1)coqvar=matrix(1,nqv,1,n);  /**< Fixed quantitative covariate */
-  if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,1,n);  /**< Time varying quantitative covariate */
-  if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,1,n);  /**< Time varying covariate (dummy and quantitative)*/
+  covar=matrix(0,NCOVMAX,firstobs,lastobs);  /**< used in readdata */
+  if(nqv>=1)coqvar=matrix(1,nqv,firstobs,lastobs);  /**< Fixed quantitative covariate */
+  if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,firstobs,lastobs);  /**< Time varying quantitative covariate */
+  if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,firstobs,lastobs);  /**< Time varying covariate (dummy and quantitative)*/
   cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/
   /* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5
      v1+v2*age+v2*v3 makes cptcovn = 3
@@ -11173,6 +11231,15 @@ int main(int argc, char *argv[])
       for(jj=1; jj <=nlstate+ndeath; jj++){
 	if(jj==i) continue;
 	j++;
+	while((c=getc(ficpar))=='#' && c!= EOF){
+	  ungetc(c,ficpar);
+	  fgets(line, MAXLINE, ficpar);
+	  numlinepar++;
+	  fputs(line,stdout);
+	  fputs(line,ficparo);
+	  fputs(line,ficlog);
+	}
+	ungetc(c,ficpar);
 	fscanf(ficpar,"%1d%1d",&i1,&j1);
 	if ((i1 != i) || (j1 != jj)){
 	  printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \
@@ -11313,16 +11380,25 @@ Please run with mle=-1 to get a correct
   
   /*  Main data
    */
-  n= lastobs;
-  num=lvector(1,n);
-  moisnais=vector(1,n);
-  annais=vector(1,n);
-  moisdc=vector(1,n);
-  andc=vector(1,n);
-  weight=vector(1,n);
-  agedc=vector(1,n);
-  cod=ivector(1,n);
-  for(i=1;i<=n;i++){
+  nobs=lastobs-firstobs+1; /* was = lastobs;*/
+  /* num=lvector(1,n); */
+  /* moisnais=vector(1,n); */
+  /* annais=vector(1,n); */
+  /* moisdc=vector(1,n); */
+  /* andc=vector(1,n); */
+  /* weight=vector(1,n); */
+  /* agedc=vector(1,n); */
+  /* cod=ivector(1,n); */
+  /* for(i=1;i<=n;i++){ */
+  num=lvector(firstobs,lastobs);
+  moisnais=vector(firstobs,lastobs);
+  annais=vector(firstobs,lastobs);
+  moisdc=vector(firstobs,lastobs);
+  andc=vector(firstobs,lastobs);
+  weight=vector(firstobs,lastobs);
+  agedc=vector(firstobs,lastobs);
+  cod=ivector(firstobs,lastobs);
+  for(i=firstobs;i<=lastobs;i++){
     num[i]=0;
     moisnais[i]=0;
     annais[i]=0;
@@ -11332,9 +11408,9 @@ Please run with mle=-1 to get a correct
     cod[i]=0;
     weight[i]=1.0; /* Equal weights, 1 by default */
   }
-  mint=matrix(1,maxwav,1,n);
-  anint=matrix(1,maxwav,1,n);
-  s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */ 
+  mint=matrix(1,maxwav,firstobs,lastobs);
+  anint=matrix(1,maxwav,firstobs,lastobs);
+  s=imatrix(1,maxwav+1,firstobs,lastobs); /* s[i][j] health state for wave i and individual j */ 
   tab=ivector(1,NCOVMAX);
   ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */
   ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */
@@ -11436,8 +11512,8 @@ Please run with mle=-1 to get a correct
 
 
   agegomp=(int)agemin;
-  free_vector(moisnais,1,n);
-  free_vector(annais,1,n);
+  free_vector(moisnais,firstobs,lastobs);
+  free_vector(annais,firstobs,lastobs);
   /* free_matrix(mint,1,maxwav,1,n);
      free_matrix(anint,1,maxwav,1,n);*/
   /* free_vector(moisdc,1,n); */
@@ -11463,8 +11539,8 @@ Please run with mle=-1 to get a correct
   concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
   /* Concatenates waves */
  
-  free_vector(moisdc,1,n);
-  free_vector(andc,1,n);
+  free_vector(moisdc,firstobs,lastobs);
+  free_vector(andc,firstobs,lastobs);
 
   /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */
   nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); 
@@ -11654,10 +11730,10 @@ Title=%s <br>Datafile=%s Firstpass=%d La
   fprintf(fichtm,"\n<li> Number of fixed quantitative variables: nqv=%d ", nqv); 
   ncurrv=i;
   for(i=ncurrv; i <=ncurrv-1+nqv; i++) fprintf(fichtm,"V%d ", i);
-  fprintf(fichtm,"\n<li> Number of time varying (wave varying) covariates: ntv=%d ", ntv);
+  fprintf(fichtm,"\n<li> Number of time varying (wave varying) dummy covariates: ntv=%d ", ntv);
   ncurrv=i;
   for(i=ncurrv; i <=ncurrv-1+ntv; i++) fprintf(fichtm,"V%d ", i);
-  fprintf(fichtm,"\n<li>Number of quantitative time varying covariates: nqtv=%d ", nqtv);
+  fprintf(fichtm,"\n<li>Number of time varying  quantitative covariates: nqtv=%d ", nqtv);
   ncurrv=i;
   for(i=ncurrv; i <=ncurrv-1+nqtv; i++) fprintf(fichtm,"V%d ", i);
   fprintf(fichtm,"\n<li>Weights column \n<br>Number of alive states: nlstate=%d <br>Number of death states (not really implemented): ndeath=%d \n<li>Number of waves: maxwav=%d \n<li>Parameter for maximization (1), using parameter values (0), for design of parameters and variance-covariance matrix: mle=%d \n<li>Does the weight column be taken into account (1), or not (0): weight=%d</ul>\n", \
@@ -11689,10 +11765,10 @@ Interval (in months) between two waves:
       for(j=1;j<=NDIM;j++)
 	ximort[i][j]=0.;
     /*     ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */
-    cens=ivector(1,n);
-    ageexmed=vector(1,n);
-    agecens=vector(1,n);
-    dcwave=ivector(1,n);
+    cens=ivector(firstobs,lastobs);
+    ageexmed=vector(firstobs,lastobs);
+    agecens=vector(firstobs,lastobs);
+    dcwave=ivector(firstobs,lastobs);
 		
     for (i=1; i<=imx; i++){
       dcwave[i]=-1;
@@ -11906,9 +11982,10 @@ Please run with mle=-1 to get a correct
     free_vector(lpop,1,AGESUP);
     free_vector(tpop,1,AGESUP);
     free_matrix(ximort,1,NDIM,1,NDIM);
-    free_ivector(cens,1,n);
-    free_vector(agecens,1,n);
-    free_ivector(dcwave,1,n);
+    free_ivector(dcwave,firstobs,lastobs);
+    free_vector(agecens,firstobs,lastobs);
+    free_vector(ageexmed,firstobs,lastobs);
+    free_ivector(cens,firstobs,lastobs);
 #ifdef GSL
 #endif
   } /* Endof if mle==-3 mortality only */
@@ -12314,8 +12391,8 @@ Please run with mle=-1 to get a correct
     /* 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);    */
-    free_lvector(num,1,n);
-    free_vector(agedc,1,n);
+    free_lvector(num,firstobs,lastobs);
+    free_vector(agedc,firstobs,lastobs);
     /*free_matrix(covar,0,NCOVMAX,1,n);*/
     /*free_matrix(covar,1,NCOVMAX,1,n);*/
     fclose(ficparo);
@@ -12572,13 +12649,13 @@ Please run with mle=-1 to get a correct
 	if(vpopbased==1)
 	  fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
 	else
-	  fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
+	  fprintf(ficrest,"the age specific forward period (stable) prevalences in each health state \n");
 	fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
 	for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
 	fprintf(ficrest,"\n");
 	/* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
-	printf("Computing age specific period (stable) prevalences in each health state \n");
-	fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
+	printf("Computing age specific forward period (stable) prevalences in each health state \n");
+	fprintf(ficlog,"Computing age specific forward period (stable) prevalences in each health state \n");
 	for(age=bage; age <=fage ;age++){
 	  prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k, nres); /*ZZ Is it the correct prevalim */
 	  if (vpopbased==1) {
@@ -12625,16 +12702,16 @@ Please run with mle=-1 to get a correct
     printf("done State-specific expectancies\n");fflush(stdout);
     fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog);
 
-    /* variance-covariance of period prevalence*/
+    /* variance-covariance of forward period prevalence*/
     varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff);
 
     
-    free_vector(weight,1,n);
+    free_vector(weight,firstobs,lastobs);
     free_imatrix(Tvard,1,NCOVMAX,1,2);
-    free_imatrix(s,1,maxwav+1,1,n);
-    free_matrix(anint,1,maxwav,1,n); 
-    free_matrix(mint,1,maxwav,1,n);
-    free_ivector(cod,1,n);
+    free_imatrix(s,1,maxwav+1,firstobs,lastobs);
+    free_matrix(anint,1,maxwav,firstobs,lastobs); 
+    free_matrix(mint,1,maxwav,firstobs,lastobs);
+    free_ivector(cod,firstobs,lastobs);
     free_ivector(tab,1,NCOVMAX);
     fclose(ficresstdeij);
     fclose(ficrescveij);
@@ -12654,10 +12731,10 @@ Please run with mle=-1 to get a correct
   free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
-  if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,1,n);
-  if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n);
-  if(nqv>=1)free_matrix(coqvar,1,nqv,1,n);
-  free_matrix(covar,0,NCOVMAX,1,n);
+  if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,firstobs,lastobs);
+  if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,firstobs,lastobs);
+  if(nqv>=1)free_matrix(coqvar,1,nqv,firstobs,lastobs);
+  free_matrix(covar,0,NCOVMAX,firstobs,lastobs);
   free_matrix(matcov,1,npar,1,npar);
   free_matrix(hess,1,npar,1,npar);
   /*free_vector(delti,1,npar);*/
@@ -12777,13 +12854,16 @@ Please run with mle=-1 to get a correct
   
   sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot);
   printf("Starting graphs with: '%s'\n",plotcmd);fflush(stdout);
+  strcpy(pplotcmd,plotcmd);
   
   if((outcmd=system(plotcmd)) != 0){
-    printf("gnuplot command might not be in your path: '%s', err=%d\n", plotcmd, outcmd);
+    printf("Error in gnuplot, command might not be in your path: '%s', err=%d\n", plotcmd, outcmd);
     printf("\n Trying if gnuplot resides on the same directory that IMaCh\n");
     sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot);
-    if((outcmd=system(plotcmd)) != 0)
+    if((outcmd=system(plotcmd)) != 0){
       printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd);
+      strcpy(plotcmd,pplotcmd);
+    }
   }
   printf(" Successful, please wait...");
   while (z[0] != 'q') {