--- imach/src/imach.c	2016/09/02 11:11:21	1.247
+++ imach/src/imach.c	2017/04/05 10:14:09	1.261
@@ -1,6 +1,50 @@
-/* $Id: imach.c,v 1.247 2016/09/02 11:11:21 brouard Exp $
+/* $Id: imach.c,v 1.261 2017/04/05 10:14:09 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.261  2017/04/05 10:14:09  brouard
+  Summary: Bug in E_ as well as in T_ fixed nres-1 vs k1-1
+
+  Revision 1.260  2017/04/04 17:46:59  brouard
+  Summary: Gnuplot indexations fixed (humm)
+
+  Revision 1.259  2017/04/04 13:01:16  brouard
+  Summary: Some errors to warnings only if date of death is unknown but status is death we could set to pi3
+
+  Revision 1.258  2017/04/03 10:17:47  brouard
+  Summary: Version 0.99r12
+
+  Some cleanings, conformed with updated documentation.
+
+  Revision 1.257  2017/03/29 16:53:30  brouard
+  Summary: Temp
+
+  Revision 1.256  2017/03/27 05:50:23  brouard
+  Summary: Temporary
+
+  Revision 1.255  2017/03/08 16:02:28  brouard
+  Summary: IMaCh version 0.99r10 bugs in gnuplot fixed
+
+  Revision 1.254  2017/03/08 07:13:00  brouard
+  Summary: Fixing data parameter line
+
+  Revision 1.253  2016/12/15 11:59:41  brouard
+  Summary: 0.99 in progress
+
+  Revision 1.252  2016/09/15 21:15:37  brouard
+  *** empty log message ***
+
+  Revision 1.251  2016/09/15 15:01:13  brouard
+  Summary: not working
+
+  Revision 1.250  2016/09/08 16:07:27  brouard
+  Summary: continue
+
+  Revision 1.249  2016/09/07 17:14:18  brouard
+  Summary: Starting values from frequencies
+
+  Revision 1.248  2016/09/07 14:10:18  brouard
+  *** empty log message ***
+
   Revision 1.247  2016/09/02 11:11:21  brouard
   *** empty log message ***
 
@@ -116,9 +160,7 @@
   Author: Nicolas Brouard
 
   Revision 1.210  2015/11/18 17:41:20  brouard
-  Summary: Start working on projected prevalences
-
-  Revision 1.209  2015/11/17 22:12:03  brouard
+  Summary: Start working on projected prevalences  Revision 1.209  2015/11/17 22:12:03  brouard
   Summary: Adding ftolpl parameter
   Author: N Brouard
 
@@ -906,7 +948,7 @@ typedef struct {
 /* #include <libintl.h> */
 /* #define _(String) gettext (String) */
 
-#define MAXLINE 1024 /* Was 256. Overflow with 312 with 2 states and 4 covariates. Should be ok */
+#define MAXLINE 2048 /* Was 256 and 1024. Overflow with 312 with 2 states and 4 covariates. Should be ok */
 
 #define GNUPLOTPROGRAM "gnuplot"
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
@@ -943,12 +985,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.247 2016/09/02 11:11:21 brouard Exp $ */
+/* $Id: imach.c,v 1.261 2017/04/05 10:14:09 brouard Exp $ */
 /* $State: Exp $ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
 char copyright[]="February 2016,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.247 $ $Date: 2016/09/02 11:11:21 $"; 
+char fullversion[]="$Revision: 1.261 $ $Date: 2017/04/05 10:14:09 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1149,6 +1191,7 @@ int *TvarsQind;
 
 #define MAXRESULTLINES 10
 int nresult=0;
+int parameterline=0; /* # of the parameter (type) line */
 int TKresult[MAXRESULTLINES];
 int Tresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
 int Tinvresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */
@@ -2381,8 +2424,8 @@ void powell(double p[], double **xi, int
 	    flatd++;
 	}
 	if(flatd >0){
-	  printf("%d flat directions\n",flatd);
-	  fprintf(ficlog,"%d flat directions\n",flatd);
+	  printf("%d flat directions: ",flatd);
+	  fprintf(ficlog,"%d flat directions :",flatd);
 	  for (j=1;j<=n;j++) { 
 	    if(flatdir[j]>0){
 	      printf("%d ",j);
@@ -4110,7 +4153,16 @@ void ludcmp(double **a, int n, int *indx
     big=0.0; 
     for (j=1;j<=n;j++) 
       if ((temp=fabs(a[i][j])) > big) big=temp; 
-    if (big == 0.0) nrerror("Singular matrix in routine ludcmp"); 
+    if (big == 0.0){
+      printf(" Singular Hessian matrix at row %d:\n",i);
+      for (j=1;j<=n;j++) {
+	printf(" a[%d][%d]=%f,",i,j,a[i][j]);
+	fprintf(ficlog," a[%d][%d]=%f,",i,j,a[i][j]);
+      }
+      fflush(ficlog);
+      fclose(ficlog);
+      nrerror("Singular matrix in routine ludcmp"); 
+    }
     vv[i]=1.0/big; 
   } 
   for (j=1;j<=n;j++) { 
@@ -4176,17 +4228,84 @@ void pstamp(FILE *fichier)
   fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart);
 }
 
+int linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb) {
+
+  /* y=a+bx regression */
+  double   sumx = 0.0;                        /* sum of x                      */
+  double   sumx2 = 0.0;                       /* sum of x**2                   */
+  double   sumxy = 0.0;                       /* sum of x * y                  */
+  double   sumy = 0.0;                        /* sum of y                      */
+  double   sumy2 = 0.0;                       /* sum of y**2                   */
+  double   sume2; /* sum of square or residuals */
+  double yhat;
+  
+  double denom=0;
+  int i;
+  int ne=*no;
+  
+  for ( i=ifi, ne=0;i<=ila;i++) {
+    if(!isfinite(x[i]) || !isfinite(y[i])){
+      /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
+      continue;
+    }
+    ne=ne+1;
+    sumx  += x[i];       
+    sumx2 += x[i]*x[i];  
+    sumxy += x[i] * y[i];
+    sumy  += y[i];      
+    sumy2 += y[i]*y[i]; 
+    denom = (ne * sumx2 - sumx*sumx);
+    /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
+  } 
+  
+  denom = (ne * sumx2 - sumx*sumx);
+  if (denom == 0) {
+    // vertical, slope m is infinity
+    *b = INFINITY;
+    *a = 0;
+    if (r) *r = 0;
+    return 1;
+  }
+  
+  *b = (ne * sumxy  -  sumx * sumy) / denom;
+  *a = (sumy * sumx2  -  sumx * sumxy) / denom;
+  if (r!=NULL) {
+    *r = (sumxy - sumx * sumy / ne) /          /* compute correlation coeff     */
+      sqrt((sumx2 - sumx*sumx/ne) *
+	   (sumy2 - sumy*sumy/ne));
+  }
+  *no=ne;
+  for ( i=ifi, ne=0;i<=ila;i++) {
+    if(!isfinite(x[i]) || !isfinite(y[i])){
+      /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
+      continue;
+    }
+    ne=ne+1;
+    yhat = y[i] - *a -*b* x[i];
+    sume2  += yhat * yhat ;       
+    
+    denom = (ne * sumx2 - sumx*sumx);
+    /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
+  } 
+  *sb = sqrt(sume2/(ne-2)/(sumx2 - sumx * sumx /ne));
+  *sa= *sb * sqrt(sumx2/ne);
+  
+  return 0; 
+}
+
 /************ Frequencies ********************/
-void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
+void  freqsummary(char fileres[], double p[], double pstart[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
 		  int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[],	\
 		  int firstpass,  int lastpass, int stepm, int weightopt, char model[])
-{  /* Some frequencies */
+{  /* Some frequencies as well as proposing some starting values */
   
-  int i, m, jk, j1, bool, z1,j, k, iv;
+  int i, m, jk, j1, bool, z1,j, nj, nl, k, iv, jj=0;
   int iind=0, iage=0;
   int mi; /* Effective wave */
   int first;
   double ***freq; /* Frequencies */
+  double *x, *y, a,b,r, sa, sb; /* for regression, y=b+m*x and r is the correlation coefficient */
+  int no;
   double *meanq;
   double **meanqt;
   double *pp, **prop, *posprop, *pospropt;
@@ -4195,7 +4314,7 @@ void  freqsummary(char fileres[], int ia
   double agebegin, ageend;
     
   pp=vector(1,nlstate);
-  prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); 
+  prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+4+AGEMARGE); 
   posprop=vector(1,nlstate); /* Counting the number of transition starting from a live state per age */ 
   pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ 
   /* prop=matrix(1,nlstate,iagemin,iagemax+3); */
@@ -4239,14 +4358,15 @@ 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 of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
   
-  freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
+  y= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE);
+  x= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE);
+  freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+4+AGEMARGE);
   j1=0;
   
   /* j=ncoveff;  /\* Only fixed dummy covariates *\/ */
   j=cptcoveff;  /* Only dummy covariates of the model */
   if (cptcovn<1) {j=1;ncodemax[1]=1;}
   
-  first=1;
   
   /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels:
      reference=low_education V1=0,V2=0
@@ -4254,277 +4374,414 @@ Title=%s <br>Datafile=%s Firstpass=%d La
      high_educ               V1=0 V2=1
      Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff 
   */
-  
-  for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives V4=0, V3=0 for example, fixed or varying covariates */
-    posproptt=0.;
-    /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
-      scanf("%d", i);*/
-    for (i=-5; i<=nlstate+ndeath; i++)  
-      for (jk=-5; jk<=nlstate+ndeath; jk++)  
+  dateintsum=0;
+  k2cpt=0;
+
+  if(cptcoveff == 0 )
+    nl=1;  /* Constant model only */
+  else
+    nl=2;
+  for (nj = 1; nj <= nl; nj++){   /* nj= 1 constant model, nl number of loops. */
+    if(nj==1)
+      j=0;  /* First pass for the constant */
+    else
+      j=cptcoveff; /* Other passes for the covariate values */
+    first=1;
+    for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */
+      posproptt=0.;
+      /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
+	scanf("%d", i);*/
+      for (i=-5; i<=nlstate+ndeath; i++)  
+	for (jk=-5; jk<=nlstate+ndeath; jk++)  
+	  for(m=iagemin; m <= iagemax+3; m++)
+	    freq[i][jk][m]=0;
+      
+      for (i=1; i<=nlstate; i++)  {
 	for(m=iagemin; m <= iagemax+3; m++)
-	  freq[i][jk][m]=0;
-    
-    for (i=1; i<=nlstate; i++)  {
-      for(m=iagemin; m <= iagemax+3; m++)
-	prop[i][m]=0;
-      posprop[i]=0;
-      pospropt[i]=0;
-    }
-    /* for (z1=1; z1<= nqfveff; z1++) {   */
-    /*   meanq[z1]+=0.; */
-    /*   for(m=1;m<=lastpass;m++){ */
-    /* 	meanqt[m][z1]=0.; */
-    /*   } */
-    /* } */
-    
-    dateintsum=0;
-    k2cpt=0;
-    /* For that combination of covariate j1, we count and print the frequencies in one pass */
-    for (iind=1; iind<=imx; iind++) { /* For each individual iind */
-      bool=1;
-      if(anyvaryingduminmodel==0){ /* If All fixed covariates */
-	if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
-	  /* for (z1=1; z1<= nqfveff; z1++) {   */
-	  /*   meanq[z1]+=coqvar[Tvar[z1]][iind];  /\* Computes mean of quantitative with selected filter *\/ */
-	  /* } */
-	  for (z1=1; z1<=cptcoveff; z1++) {  
-	    /* if(Tvaraff[z1] ==-20){ */
-	    /* 	 /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
-	    /* }else  if(Tvaraff[z1] ==-10){ */
-	    /* 	 /\* sumnew+=coqvar[z1][iind]; *\/ */
-	    /* }else  */
-	    if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
-	      /* Tests if this individual iind responded to j1 (V4=1 V3=0) */
-	      bool=0;
-	      /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", 
-		 bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
-		 j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
-	      /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
-	    } /* Onlyf fixed */
-	  } /* end z1 */
-	} /* cptcovn > 0 */
-      } /* end any */
-      if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */
-	/* for(m=firstpass; m<=lastpass; m++){ */
-	for(mi=1; mi<wav[iind];mi++){ /* For that wave */
-	  m=mw[mi][iind];
-	  if(anyvaryingduminmodel==1){ /* Some are varying covariates */
-	    for (z1=1; z1<=cptcoveff; z1++) {
-	      if( Fixed[Tmodelind[z1]]==1){
-		iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
-		if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality */
-		  bool=0;
-	      }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */
-		if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) {
-		  bool=0;
+	  prop[i][m]=0;
+	posprop[i]=0;
+	pospropt[i]=0;
+      }
+      /* for (z1=1; z1<= nqfveff; z1++) {   */
+      /*   meanq[z1]+=0.; */
+      /*   for(m=1;m<=lastpass;m++){ */
+      /* 	meanqt[m][z1]=0.; */
+      /*   } */
+      /* } */
+      
+      /* dateintsum=0; */
+      /* k2cpt=0; */
+      
+      /* For that combination of covariate j1, we count and print the frequencies in one pass */
+      for (iind=1; iind<=imx; iind++) { /* For each individual iind */
+	bool=1;
+	if(j !=0){
+	  if(anyvaryingduminmodel==0){ /* If All fixed covariates */
+	    if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
+	      /* for (z1=1; z1<= nqfveff; z1++) {   */
+	      /*   meanq[z1]+=coqvar[Tvar[z1]][iind];  /\* Computes mean of quantitative with selected filter *\/ */
+	      /* } */
+	      for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */
+		/* if(Tvaraff[z1] ==-20){ */
+		/* 	 /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
+		/* }else  if(Tvaraff[z1] ==-10){ */
+		/* 	 /\* sumnew+=coqvar[z1][iind]; *\/ */
+		/* }else  */
+		if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ /* for combination j1 of covariates */
+		  /* Tests if this individual iind responded to combination j1 (V4=1 V3=0) */
+		  bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */
+		  /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", 
+		     bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
+		     j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
+		  /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
+		} /* Onlyf fixed */
+	      } /* end z1 */
+	    } /* cptcovn > 0 */
+	  } /* end any */
+	}/* end j==0 */
+	if (bool==1){ /* We selected an individual iind satisfying combination j1 or all fixed */
+	  /* for(m=firstpass; m<=lastpass; m++){ */
+	  for(mi=1; mi<wav[iind];mi++){ /* For that wave */
+	    m=mw[mi][iind];
+	    if(j!=0){
+	      if(anyvaryingduminmodel==1){ /* Some are varying covariates */
+		for (z1=1; z1<=cptcoveff; z1++) {
+		  if( Fixed[Tmodelind[z1]]==1){
+		    iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
+		    if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) /* iv=1 to ntv, right modality. If covariate's 
+										      value is -1, we don't select. It differs from the 
+										      constant and age model which counts them. */
+		      bool=0; /* not selected */
+		  }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */
+		    if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) {
+		      bool=0;
+		    }
+		  }
 		}
+	      }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop  */
+	    } /* end j==0 */
+	    /* bool =0 we keep that guy which corresponds to the combination of dummy values */
+	    if(bool==1){
+	      /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]
+		 and mw[mi+1][iind]. dh depends on stepm. */
+	      agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/
+	      ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */
+	      if(m >=firstpass && m <=lastpass){
+		k2=anint[m][iind]+(mint[m][iind]/12.);
+		/*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
+		if(agev[m][iind]==0) agev[m][iind]=iagemax+1;  /* All ages equal to 0 are in iagemax+1 */
+		if(agev[m][iind]==1) agev[m][iind]=iagemax+2;  /* All ages equal to 1 are in iagemax+2 */
+		if (s[m][iind]>0 && s[m][iind]<=nlstate)  /* If status at wave m is known and a live state */
+		  prop[s[m][iind]][(int)agev[m][iind]] += weight[iind];  /* At age of beginning of transition, where status is known */
+		if (m<lastpass) {
+		  /* if(s[m][iind]==4 && s[m+1][iind]==4) */
+		  /*   printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */
+		  if(s[m][iind]==-1)
+		    printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));
+		  freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
+		  /* if((int)agev[m][iind] == 55) */
+		  /*   printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */
+		  /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */
+		  freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
+		}
+	      } /* end if between passes */  
+	      if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99) && (j==0)) {
+		dateintsum=dateintsum+k2; /* on all covariates ?*/
+		k2cpt++;
+		/* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
 	      }
-	    }
-	  }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop  */
-	  /* bool =0 we keep that guy which corresponds to the combination of dummy values */
-	  if(bool==1){
-	    /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]
-	       and mw[mi+1][iind]. dh depends on stepm. */
-	    agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/
-	    ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */
-	    if(m >=firstpass && m <=lastpass){
-	      k2=anint[m][iind]+(mint[m][iind]/12.);
-	      /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
-	      if(agev[m][iind]==0) agev[m][iind]=iagemax+1;  /* All ages equal to 0 are in iagemax+1 */
-	      if(agev[m][iind]==1) agev[m][iind]=iagemax+2;  /* All ages equal to 1 are in iagemax+2 */
-	      if (s[m][iind]>0 && s[m][iind]<=nlstate)  /* If status at wave m is known and a live state */
-		prop[s[m][iind]][(int)agev[m][iind]] += weight[iind];  /* At age of beginning of transition, where status is known */
-	      if (m<lastpass) {
-		/* if(s[m][iind]==4 && s[m+1][iind]==4) */
-		/*   printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */
-		if(s[m][iind]==-1)
-		  printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));
-		freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
-		/* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */
-		freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
-	      }
-	    } /* end if between passes */  
-	    if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) {
-	      dateintsum=dateintsum+k2;
-	      k2cpt++;
-	      /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
-	    }
-	  } /* end bool 2 */
-	} /* end m */
-      } /* end bool */
-    } /* end iind = 1 to imx */
-    /* prop[s][age] is feeded for any initial and valid live state as well as
-       freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */
-    
-    
-    /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
-    pstamp(ficresp);
-    if  (cptcoveff>0){
-      fprintf(ficresp, "\n#********** Variable "); 
-      fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); 
-      fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); 
-      fprintf(ficlog, "\n#********** Variable "); 
-      for (z1=1; z1<=cptcoveff; z1++){
-	if(DummyV[z1]){
-	  fprintf(ficresp, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-	  fprintf(ficresphtm, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-	  fprintf(ficresphtmfr, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-	  fprintf(ficlog, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-	}else{
-	  fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-	  fprintf(ficresphtm, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-	  fprintf(ficresphtmfr, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-	  fprintf(ficlog, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-	}
-      }
-      fprintf(ficresp, "**********\n#");
-      fprintf(ficresphtm, "**********</h3>\n");
-      fprintf(ficresphtmfr, "**********</h3>\n");
-      fprintf(ficlog, "**********\n");
-    }
-    fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
-    for(i=1; i<=nlstate;i++) {
-      fprintf(ficresp, " Age Prev(%d)  N(%d)  N  ",i,i);
-      fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);
-    }
-    fprintf(ficresp, "\n");
-    fprintf(ficresphtm, "\n");
-    
-    /* Header of frequency table by age */
-    fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");
-    fprintf(ficresphtmfr,"<th>Age</th> ");
-    for(jk=-1; jk <=nlstate+ndeath; jk++){
-      for(m=-1; m <=nlstate+ndeath; m++){
-	if(jk!=0 && m!=0)
-	  fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);
-      }
-    }
-    fprintf(ficresphtmfr, "\n");
-    
-    /* For each age */
-    for(iage=iagemin; iage <= iagemax+3; iage++){
-      fprintf(ficresphtm,"<tr>");
-      if(iage==iagemax+1){
-	fprintf(ficlog,"1");
-	fprintf(ficresphtmfr,"<tr><th>0</th> ");
-      }else if(iage==iagemax+2){
-	fprintf(ficlog,"0");
-	fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
-      }else if(iage==iagemax+3){
-	fprintf(ficlog,"Total");
-	fprintf(ficresphtmfr,"<tr><th>Total</th> ");
-      }else{
-	if(first==1){
-	  first=0;
-	  printf("See log file for details...\n");
+	    }else{
+	      bool=1;
+	    }/* end bool 2 */
+	  } /* end m */
+	} /* end bool */
+      } /* end iind = 1 to imx */
+      /* prop[s][age] is feeded for any initial and valid live state as well as
+	 freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */
+      
+      
+      /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
+      pstamp(ficresp);
+      if  (cptcoveff>0 && j!=0){
+	printf( "\n#********** Variable "); 
+	fprintf(ficresp, "\n#********** Variable "); 
+	fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); 
+	fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); 
+	fprintf(ficlog, "\n#********** Variable "); 
+	for (z1=1; z1<=cptcoveff; z1++){
+	  if(!FixedV[Tvaraff[z1]]){
+	    printf( "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	    fprintf(ficresp, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	    fprintf(ficresphtm, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	    fprintf(ficresphtmfr, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	    fprintf(ficlog, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	  }else{
+	    printf( "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	    fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	    fprintf(ficresphtm, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	    fprintf(ficresphtmfr, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	    fprintf(ficlog, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	  }
+	}
+	printf( "**********\n#");
+	fprintf(ficresp, "**********\n#");
+	fprintf(ficresphtm, "**********</h3>\n");
+	fprintf(ficresphtmfr, "**********</h3>\n");
+	fprintf(ficlog, "**********\n");
+      }
+      fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
+      for(i=1; i<=nlstate;i++) {
+	fprintf(ficresp, " Age Prev(%d)  N(%d)  N  ",i,i);
+	fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);
+      }
+      fprintf(ficresp, "\n");
+      fprintf(ficresphtm, "\n");
+      
+      /* Header of frequency table by age */
+      fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");
+      fprintf(ficresphtmfr,"<th>Age</th> ");
+      for(jk=-1; jk <=nlstate+ndeath; jk++){
+	for(m=-1; m <=nlstate+ndeath; m++){
+	  if(jk!=0 && m!=0)
+	    fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);
 	}
-	fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage);
-	fprintf(ficlog,"Age %d", iage);
       }
-      for(jk=1; jk <=nlstate ; jk++){
-	for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
-	  pp[jk] += freq[jk][m][iage]; 
-      }
-      for(jk=1; jk <=nlstate ; jk++){
-	for(m=-1, pos=0; m <=0 ; m++)
-	  pos += freq[jk][m][iage];
-	if(pp[jk]>=1.e-10){
+      fprintf(ficresphtmfr, "\n");
+    
+      /* For each age */
+      for(iage=iagemin; iage <= iagemax+3; iage++){
+	fprintf(ficresphtm,"<tr>");
+	if(iage==iagemax+1){
+	  fprintf(ficlog,"1");
+	  fprintf(ficresphtmfr,"<tr><th>0</th> ");
+	}else if(iage==iagemax+2){
+	  fprintf(ficlog,"0");
+	  fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
+	}else if(iage==iagemax+3){
+	  fprintf(ficlog,"Total");
+	  fprintf(ficresphtmfr,"<tr><th>Total</th> ");
+	}else{
 	  if(first==1){
-	    printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
+	    first=0;
+	    printf("See log file for details...\n");
+	  }
+	  fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage);
+	  fprintf(ficlog,"Age %d", iage);
+	}
+	for(jk=1; jk <=nlstate ; jk++){
+	  for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
+	    pp[jk] += freq[jk][m][iage]; 
+	}
+	for(jk=1; jk <=nlstate ; jk++){
+	  for(m=-1, pos=0; m <=0 ; m++)
+	    pos += freq[jk][m][iage];
+	  if(pp[jk]>=1.e-10){
+	    if(first==1){
+	      printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
+	    }
+	    fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
+	  }else{
+	    if(first==1)
+	      printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
+	    fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
 	  }
-	  fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
-	}else{
-	  if(first==1)
-	    printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
-	  fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
 	}
-      }
-      
-      for(jk=1; jk <=nlstate ; jk++){ 
-	/* posprop[jk]=0; */
-	for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */
-	  pp[jk] += freq[jk][m][iage];
-      }	/* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */
       
-      for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){
-	pos += pp[jk]; /* pos is the total number of transitions until this age */
-	posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state
+	for(jk=1; jk <=nlstate ; jk++){ 
+	  /* posprop[jk]=0; */
+	  for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */
+	    pp[jk] += freq[jk][m][iage];
+	}	/* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */
+      
+	for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){
+	  pos += pp[jk]; /* pos is the total number of transitions until this age */
+	  posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state
+					    from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
+	  pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state
 					  from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
-	pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state
-					from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
-      }
-      for(jk=1; jk <=nlstate ; jk++){
-	if(pos>=1.e-5){
-	  if(first==1)
-	    printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
-	  fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
-	}else{
-	  if(first==1)
-	    printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
-	  fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
 	}
-	if( iage <= iagemax){
+	for(jk=1; jk <=nlstate ; jk++){
 	  if(pos>=1.e-5){
-	    fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
-	    fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
-	    /*probs[iage][jk][j1]= pp[jk]/pos;*/
-	    /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/
+	    if(first==1)
+	      printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
+	    fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
+	  }else{
+	    if(first==1)
+	      printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
+	    fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
+	  }
+	  if( iage <= iagemax){
+	    if(pos>=1.e-5){
+	      fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
+	      fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
+	      /*probs[iage][jk][j1]= pp[jk]/pos;*/
+	      /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/
+	    }
+	    else{
+	      fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);
+	      fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);
+	    }
 	  }
-	  else{
-	    fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);
-	    fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);
+	  pospropt[jk] +=posprop[jk];
+	} /* end loop jk */
+	/* pospropt=0.; */
+	for(jk=-1; jk <=nlstate+ndeath; jk++){
+	  for(m=-1; m <=nlstate+ndeath; m++){
+	    if(freq[jk][m][iage] !=0 ) { /* minimizing output */
+	      if(first==1){
+		printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
+	      }
+	      /* printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); */
+	      fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
+	    }
+	    if(jk!=0 && m!=0)
+	      fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]);
 	  }
+	} /* end loop jk */
+	posproptt=0.; 
+	for(jk=1; jk <=nlstate; jk++){
+	  posproptt += pospropt[jk];
+	}
+	fprintf(ficresphtmfr,"</tr>\n ");
+	if(iage <= iagemax){
+	  fprintf(ficresp,"\n");
+	  fprintf(ficresphtm,"</tr>\n");
 	}
-	pospropt[jk] +=posprop[jk];
-      } /* end loop jk */
-      /* pospropt=0.; */
-      for(jk=-1; jk <=nlstate+ndeath; jk++){
-	for(m=-1; m <=nlstate+ndeath; m++){
-	  if(freq[jk][m][iage] !=0 ) { /* minimizing output */
-	    if(first==1){
-	      printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
+	if(first==1)
+	  printf("Others in log...\n");
+	fprintf(ficlog,"\n");
+      } /* end loop age iage */
+      fprintf(ficresphtm,"<tr><th>Tot</th>");
+      for(jk=1; jk <=nlstate ; jk++){
+	if(posproptt < 1.e-5){
+	  fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt);	
+	}else{
+	  fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);	
+	}
+      }
+      fprintf(ficresphtm,"</tr>\n");
+      fprintf(ficresphtm,"</table>\n");
+      fprintf(ficresphtmfr,"</table>\n");
+      if(posproptt < 1.e-5){
+	fprintf(ficresphtm,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);
+	fprintf(ficresphtmfr,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);
+	fprintf(ficlog,"#  This combination (%d) is not valid and no result will be produced\n",j1);
+	printf("#  This combination (%d) is not valid and no result will be produced\n",j1);
+	invalidvarcomb[j1]=1;
+      }else{
+	fprintf(ficresphtm,"\n <p> This combination (%d) is valid and result will be produced.</p>",j1);
+	invalidvarcomb[j1]=0;
+      }
+      fprintf(ficresphtmfr,"</table>\n");
+      fprintf(ficlog,"\n");
+      if(j!=0){
+	printf("#Freqsummary: Starting values for combination j1=%d:\n", j1);
+	for(i=1,jk=1; i <=nlstate; i++){
+	  for(k=1; k <=(nlstate+ndeath); k++){
+	    if (k != i) {
+	      for(jj=1; jj <=ncovmodel; jj++){ /* For counting jk */
+		if(jj==1){  /* Constant case (in fact cste + age) */
+		  if(j1==1){ /* All dummy covariates to zero */
+		    freq[i][k][iagemax+4]=freq[i][k][iagemax+3]; /* Stores case 0 0 0 */
+		    freq[i][i][iagemax+4]=freq[i][i][iagemax+3]; /* Stores case 0 0 0 */
+		    printf("%d%d ",i,k);
+		    fprintf(ficlog,"%d%d ",i,k);
+		    printf("%12.7f ln(%.0f/%.0f)= %f, OR=%f sd=%f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]),freq[i][k][iagemax+3]/freq[i][i][iagemax+3], sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]));
+		    fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
+		    pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]);
+		  }
+		}else if((j1==1) && (jj==2 || nagesqr==1)){ /* age or age*age parameter without covariate V4*age (to be done later) */
+		  for(iage=iagemin; iage <= iagemax+3; iage++){
+		    x[iage]= (double)iage;
+		    y[iage]= log(freq[i][k][iage]/freq[i][i][iage]);
+		    /* printf("i=%d, k=%d, jk=%d, j1=%d, jj=%d, y[%d]=%f\n",i,k,jk,j1,jj, iage, y[iage]); */
+		  }
+		  linreg(iagemin,iagemax,&no,x,y,&a,&b,&r, &sa, &sb ); /* y= a+b*x with standard errors */
+		  pstart[jk]=b;
+		  pstart[jk-1]=a;
+		}else if( j1!=1 && (j1==2 || (log(j1-1.)/log(2.)-(int)(log(j1-1.)/log(2.))) <0.010) && ( TvarsDind[(int)(log(j1-1.)/log(2.))+1]+2+nagesqr == jj)  && Dummy[jj-2-nagesqr]==0){ /* We want only if the position, jj, in model corresponds to unique covariate equal to 1 in j1 combination */ 
+		  printf("j1=%d, jj=%d, (int)(log(j1-1.)/log(2.))+1=%d, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(int)(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]);
+		  printf("j1=%d, jj=%d, (log(j1-1.)/log(2.))+1=%f, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]);
+		  pstart[jk]= log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]));
+		  printf("%d%d ",i,k);
+		  fprintf(ficlog,"%d%d ",i,k);
+		  printf("jk=%d,i=%d,k=%d,p[%d]=%12.7f ln((%.0f/%.0f)/(%.0f/%.0f))= %f, OR=%f sd=%f \n",jk,i,k,jk,p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3],freq[i][k][iagemax+4],freq[i][i][iagemax+4], log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])),(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]), sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]+1/freq[i][k][iagemax+4]+1/freq[i][i][iagemax+4]));
+		}else{ /* Other cases, like quantitative fixed or varying covariates */
+		  ;
+		}
+		/* printf("%12.7f )", param[i][jj][k]); */
+		/* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */
+		jk++; 
+	      } /* end jj */
+	    } /* end k!= i */
+	  } /* end k */
+	} /* end i, jk */
+      } /* end j !=0 */
+    } /* end selected combination of covariate j1 */
+    if(j==0){ /* We can estimate starting values from the occurences in each case */
+      printf("#Freqsummary: Starting values for the constants:\n");
+      fprintf(ficlog,"\n");
+      for(i=1,jk=1; i <=nlstate; i++){
+	for(k=1; k <=(nlstate+ndeath); k++){
+	  if (k != i) {
+	    printf("%d%d ",i,k);
+	    fprintf(ficlog,"%d%d ",i,k);
+	    for(jj=1; jj <=ncovmodel; jj++){
+	      pstart[jk]=p[jk]; /* Setting pstart to p values by default */
+	      if(jj==1){ /* Age has to be done */
+		pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]);
+		printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
+		fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
+	      }
+	      /* printf("%12.7f )", param[i][jj][k]); */
+	      /* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */
+	      jk++; 
 	    }
-	    fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
+	    printf("\n");
+	    fprintf(ficlog,"\n");
 	  }
-	  if(jk!=0 && m!=0)
-	    fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]);
 	}
-      } /* end loop jk */
-      posproptt=0.; 
-      for(jk=1; jk <=nlstate; jk++){
-	posproptt += pospropt[jk];
-      }
-      fprintf(ficresphtmfr,"</tr>\n ");
-      if(iage <= iagemax){
-	fprintf(ficresp,"\n");
-	fprintf(ficresphtm,"</tr>\n");
       }
-      if(first==1)
-	printf("Others in log...\n");
+      printf("#Freqsummary\n");
       fprintf(ficlog,"\n");
-    } /* end loop age iage */
-    fprintf(ficresphtm,"<tr><th>Tot</th>");
-    for(jk=1; jk <=nlstate ; jk++){
-      if(posproptt < 1.e-5){
-	fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt);	
-      }else{
-	fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);	
+      for(jk=-1; jk <=nlstate+ndeath; jk++){
+	for(m=-1; m <=nlstate+ndeath; m++){
+	  /* param[i]|j][k]= freq[jk][m][iagemax+3] */
+	  printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]);
+	  fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]);
+	  /* if(freq[jk][m][iage] !=0 ) { /\* minimizing output *\/ */
+	  /*   printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */
+	  /*   fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */
+	  /* } */
+	}
+      } /* end loop jk */
+      
+      printf("\n");
+      fprintf(ficlog,"\n");
+    } /* end j=0 */
+  } /* end j */
+
+  if(mle == -2){  /* We want to use these values as starting values */
+    for(i=1, jk=1; i <=nlstate; i++){
+      for(j=1; j <=nlstate+ndeath; j++){
+	if(j!=i){
+	  /*ca[0]= k+'a'-1;ca[1]='\0';*/
+	  printf("%1d%1d",i,j);
+	  fprintf(ficparo,"%1d%1d",i,j);
+	  for(k=1; k<=ncovmodel;k++){
+	    /* 	  printf(" %lf",param[i][j][k]); */
+	    /* 	  fprintf(ficparo," %lf",param[i][j][k]); */
+	    p[jk]=pstart[jk];
+	    printf(" %f ",pstart[jk]);
+	    fprintf(ficparo," %f ",pstart[jk]);
+	    jk++;
+	  }
+	  printf("\n");
+	  fprintf(ficparo,"\n");
+	}
       }
     }
-    fprintf(ficresphtm,"</tr>\n");
-    fprintf(ficresphtm,"</table>\n");
-    fprintf(ficresphtmfr,"</table>\n");
-    if(posproptt < 1.e-5){
-      fprintf(ficresphtm,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);
-      fprintf(ficresphtmfr,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);
-      fprintf(ficres,"\n  This combination (%d) is not valid and no result will be produced\n\n",j1);
-      invalidvarcomb[j1]=1;
-    }else{
-      fprintf(ficresphtm,"\n <p> This combination (%d) is valid and result will be produced.</p>",j1);
-      invalidvarcomb[j1]=0;
-    }
-    fprintf(ficresphtmfr,"</table>\n");
-  } /* end selected combination of covariate j1 */
+  } /* end mle=-2 */
   dateintmean=dateintsum/k2cpt; 
   
   fclose(ficresp);
@@ -4532,10 +4789,12 @@ Title=%s <br>Datafile=%s Firstpass=%d La
   fclose(ficresphtmfr);
   free_vector(meanq,1,nqfveff);
   free_matrix(meanqt,1,lastpass,1,nqtveff);
-  free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+3+AGEMARGE);
+  free_vector(x, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
+  free_vector(y, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
+  free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
   free_vector(pospropt,1,nlstate);
   free_vector(posprop,1,nlstate);
-  free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE);
+  free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+4+AGEMARGE);
   free_vector(pp,1,nlstate);
   /* End of freqsummary */
 }
@@ -4562,7 +4821,7 @@ void prevalence(double ***probs, double
   iagemin= (int) agemin;
   iagemax= (int) agemax;
   /*pp=vector(1,nlstate);*/
-  prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); 
+  prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+4+AGEMARGE); 
   /*  freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/
   j1=0;
   
@@ -4572,7 +4831,7 @@ void prevalence(double ***probs, double
   first=1;
   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+3+AGEMARGE; iage++)
+      for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++)
 	prop[i][iage]=0.0;
     printf("Prevalence combination of varying and fixed dummies %d\n",j1);
     /* fprintf(ficlog," V%d=%d ",Tvaraff[j1],nbcode[Tvaraff[j1]][codtabm(k,j1)]); */
@@ -4603,7 +4862,7 @@ void prevalence(double ***probs, double
 	    if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
 	      if(agev[m][i]==0) agev[m][i]=iagemax+1;
 	      if(agev[m][i]==1) agev[m][i]=iagemax+2;
-	      if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+3+AGEMARGE){
+	      if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+4+AGEMARGE){
 		printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d  m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m); 
 		exit(1);
 	      }
@@ -4640,7 +4899,7 @@ 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-AGEMARGE,iagemax+3+AGEMARGE);
+  free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+4+AGEMARGE);
 }  /* End of prevalence */
 
 /************* Waves Concatenation ***************/
@@ -4719,13 +4978,9 @@ void  concatwav(int wav[], int **dh, int
       /* if(mi==0)  never been interviewed correctly before death */
       /* Only death is a correct wave */
       mw[mi][i]=m;
-    }
+    } /* else not in a death state */
 #ifndef DISPATCHINGKNOWNDEATHAFTERLASTWAVE
-    else if ((int) andc[i] != 9999) { /* Status is negative. A death occured after lastpass, we can't take it into account because of potential bias */
-      /* m++; */
-      /* mi++; */
-      /* s[m][i]=nlstate+1;  /\* We are setting the status to the last of non live state *\/ */
-      /* mw[mi][i]=m; */
+    else if ((int) andc[i] != 9999) {  /* Date of death is known */
       if ((int)anint[m][i]!= 9999) { /* date of last interview is known */
 	if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* death occured before last wave and status should have been death instead of -1 */
 	  nbwarn++;
@@ -4738,12 +4993,12 @@ void  concatwav(int wav[], int **dh, int
 	}else{ /* Death occured afer last wave potential bias */
 	  nberr++;
 	  if(firstwo==0){
-	    printf("Error! Death for individual %ld line=%d occurred at %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 );
+	    printf("Error! Death for individual %ld line=%d occurred at %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. Please add a new fictive wave at the date of last vital status scan, with a dead status or alive but unknown state status (-1). See documentation\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
 	    firstwo=1;
 	  }
-	  fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %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,"Error! Death for individual %ld line=%d occurred at %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. Please add a new fictive wave at the date of last vital status scan, with a dead status or alive but unknown state status (-1). See documentation\n\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
 	}
-      }else{ /* end date of interview is known */
+      }else{ /* if date of interview is unknown */
 	/* 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 );
@@ -4774,7 +5029,7 @@ void  concatwav(int wav[], int **dh, int
       if (stepm <=0)
 	dh[mi][i]=1;
       else{
-	if (s[mw[mi+1][i]][i] > nlstate) { /* A death */
+	if (s[mw[mi+1][i]][i] > nlstate) { /* A death, but what if date is unknown? */
 	  if (agedc[i] < 2*AGESUP) {
 	    j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); 
 	    if(j==0) j=1;  /* Survives at least one month after exam */
@@ -6136,7 +6391,7 @@ 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 backcast, int estepm , \
+		  int popforecast, int mobilav, int prevfcast, int mobilavproj, 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, k4, nres;
@@ -6183,7 +6438,7 @@ void printinghtml(char fileresu[], char
 
    for(nres=1; nres <= nresult; nres++) /* For each resultline */
    for(k1=1; k1<=m;k1++){ /* For each combination of covariate */
-     if(TKresult[nres]!= k1)
+     if(m != 1 && TKresult[nres]!= k1)
        continue;
 
      /* for(i1=1; i1<=ncodemax[k1];i1++){ */
@@ -6210,7 +6465,7 @@ void printinghtml(char fileresu[], char
        }
      }
      /* aij, bij */
-     fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1-%d.svg\">%s_%d-1-%d.svg</a><br> \
+     fprintf(fichtm,"<br>- Logit model (yours is: logit(pij)=log(pij/pii)= aij+ bij age+%s) as a function of age: <a href=\"%s_%d-1-%d.svg\">%s_%d-1-%d.svg</a><br> \
 <img src=\"%s_%d-1-%d.svg\">",model,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
      /* Pij */
      fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2-%d.svg\">%s_%d-2-%d.svg</a><br> \
@@ -6234,21 +6489,21 @@ 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-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+       fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d some years earlier, knowing that we will be 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,"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 */
        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-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+	 fprintf(fichtm,"<br>\n- Convergence to mixed (stable) back prevalence in state %d. Or probability to be in state %d at a younger age, knowing that we will be 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 */
        for(cpt=1; cpt<=nlstate;cpt++){
-	 fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
-<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
+	 fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
        }
      }
 	 
@@ -6315,7 +6570,7 @@ See page 'Matrix of variance-covariance
 
    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
    for(k1=1; k1<=m;k1++){
-     if(TKresult[nres]!= k1)
+     if(m != 1 && TKresult[nres]!= k1)
        continue;
      /* for(i1=1; i1<=ncodemax[k1];i1++){ */
      jj1++;
@@ -6336,9 +6591,9 @@ See page 'Matrix of variance-covariance
        }
      }
      for(cpt=1; cpt<=nlstate;cpt++) {
-       fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \
+       fprintf(fichtm,"\n<br>- Observed (cross-sectional with mov_average=%d) and period (incidence based) \
 prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d-%d.svg\"> %s_%d-%d-%d.svg</a>\n <br>\
-<img src=\"%s_%d-%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres);  
+<img src=\"%s_%d-%d-%d.svg\">",mobilav,cpt,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres);  
      }
      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) \
@@ -6411,7 +6666,7 @@ void printinggnuplot(char fileresu[], ch
     for (k1=1; k1<= m ; k1 ++){ /* For each valid combination of covariate */
       for(nres=1; nres <= nresult; nres++){ /* For each resultline */
 	/* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
-	if(TKresult[nres]!= k1)
+	if(m != 1 && TKresult[nres]!= k1)
 	  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); */
@@ -6433,24 +6688,26 @@ void printinggnuplot(char fileresu[], ch
 	/* printf("\n#\n"); */
 	fprintf(ficgp,"\n#\n");
 	if(invalidvarcomb[k1]){
+          /*k1=k1-1;*/ /* To be checked */
 	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
 	  continue;
 	}
       
 	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
 	fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);
-	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres);
-      
+	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres);
+	/* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */
+      /* k1-1 error should be nres-1*/
 	for (i=1; i<= nlstate ; i ++) {
 	  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_"),k1-1,k1-1,nres);
+	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);
 	for (i=1; i<= nlstate ; i ++) {
 	  if (i==cpt) 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==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); 
+	fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%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)");
@@ -6493,7 +6750,7 @@ void printinggnuplot(char fileresu[], ch
   /*2 eme*/
   for (k1=1; k1<= m ; k1 ++){  
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
 	continue;
       fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
@@ -6523,20 +6780,20 @@ void printinggnuplot(char fileresu[], ch
 	  fprintf(ficgp,"\nreplot ");
 	for (i=1; i<= nlstate+1 ; 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_"),nres-1,nres-1, vpopbased);
 	  for (j=1; j<= nlstate+1 ; j ++) {
 	    if (j==i) fprintf(ficgp," %%lf (%%lf)");
 	    else fprintf(ficgp," %%*lf (%%*lf)");
 	  }   
 	  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);
-	  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_"),nres-1,nres-1,vpopbased);
 	  for (j=1; j<= nlstate+1 ; j ++) {
 	    if (j==i) fprintf(ficgp," %%lf (%%lf)");
 	    else fprintf(ficgp," %%*lf (%%*lf)");
 	  }   
 	  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_"),nres-1,nres-1,vpopbased);
 	  for (j=1; j<= nlstate+1 ; j ++) {
 	    if (j==i) fprintf(ficgp," %%lf (%%lf)");
 	    else fprintf(ficgp," %%*lf (%%*lf)");
@@ -6553,11 +6810,11 @@ void printinggnuplot(char fileresu[], ch
   /*3eme*/
   for (k1=1; k1<= m ; k1 ++){
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
 	continue;
 
       for (cpt=1; cpt<= nlstate ; cpt ++) {
-	fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files:  combination=%d state=%d",k1, cpt);
+	fprintf(ficgp,"\n\n# 3d: Life expectancy with EXP_ files:  combination=%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 */
@@ -6579,7 +6836,7 @@ void printinggnuplot(char fileresu[], ch
 	k=2+(nlstate+1)*(cpt-1);
 	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres);
 	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);
+plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),nres-1,nres-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) ");
 	  fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
@@ -6589,11 +6846,11 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 				
 	*/
 	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_"),nres-1,nres-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.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
+	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),nres-1,nres-1,k+nlstate,cpt);
       }
     } /* end nres */
   } /* end kl 3eme */
@@ -6602,7 +6859,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
   /* Survival functions (period) from state i in state j by initial state i */
   for (k1=1; k1<=m; k1++){    /* For each covariate and each value */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
 	continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state cpt*/
 	fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
@@ -6648,7 +6905,7 @@ set ter svg size 640, 480\nunset log y\n
   /* Survival functions (period) from state i in state j by final state j */
   for (k1=1; k1<= m ; k1++){ /* For each covariate combination if any */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
 	continue;
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
 	fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
@@ -6702,9 +6959,9 @@ set ter svg size 640, 480\nunset log y\n
   /* CV preval stable (period) 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(TKresult[nres]!= k1)
+    if(m != 1 && TKresult[nres]!= k1)
       continue;
-    for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
+    for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state of arrival */
       
       fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
@@ -6728,12 +6985,12 @@ set ter svg size 640, 480\nunset log y\n
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
 set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
       k=3; /* Offset */
-      for (i=1; i<= nlstate ; i ++){
+      for (i=1; i<= nlstate ; i ++){ /* State of origin */
 	if(i==1)
 	  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
 	else
 	  fprintf(ficgp,", '' ");
-	l=(nlstate+ndeath)*(i-1)+1;
+	l=(nlstate+ndeath)*(i-1)+1; /* 1, 1+ nlstate+ndeath, 1+2*(nlstate+ndeath) */
 	fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
 	for (j=2; j<= nlstate ; j ++)
 	  fprintf(ficgp,"+$%d",k+l+j-1);
@@ -6749,10 +7006,10 @@ set ter svg size 640, 480\nunset log y\n
     /* CV back preval stable (period) 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(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
 	continue;
-      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
-	fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
+      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life ending state */
+	fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): '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 */
@@ -6774,16 +7031,16 @@ set ter svg size 640, 480\nunset log y\n
 	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
 set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
 	k=3; /* Offset */
-	for (i=1; i<= nlstate ; i ++){
+	for (i=1; i<= nlstate ; i ++){ /* State of origin */
 	  if(i==1)
 	    fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));
 	  else
 	    fprintf(ficgp,", '' ");
 	  /* l=(nlstate+ndeath)*(i-1)+1; */
-	  l=(nlstate+ndeath)*(cpt-1)+1;
+	  l=(nlstate+ndeath)*(cpt-1)+1; /* fixed for i; cpt=1 1, cpt=2 1+ nlstate+ndeath, 1+2*(nlstate+ndeath) */
 	  /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */
 	  /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */
-	  fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+(cpt-1)+i-1); /* a vérifier */
+	  fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+i-1); /* To be verified */
 	  /* for (j=2; j<= nlstate ; j ++) */
 	  /* 	fprintf(ficgp,"+$%d",k+l+j-1); */
 	  /* 	/\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */
@@ -6800,7 +7057,7 @@ set ter svg size 640, 480\nunset log y\n
     
     for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k1)
+      if(m != 1 && TKresult[nres]!= k1)
 	continue;
       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);
@@ -6931,7 +7188,7 @@ set ter svg size 640, 480\nunset log y\n
     fprintf(ficgp,"#   jk=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */
     for(jk=1; jk <=m; jk++)  /* For each combination of covariate */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= jk)
+      if(m != 1 && TKresult[nres]!= jk)
 	continue;
       fprintf(ficgp,"# Combination of dummy  jk=%d and ",jk);
       for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
@@ -7293,7 +7550,7 @@ set ter svg size 640, 480\nunset log y\n
 /* 	      if (h==(int)(YEARM*yearp)){ */
   for(nres=1; nres <= nresult; nres++) /* For each resultline */
   for(k=1; k<=i1;k++){
-    if(TKresult[nres]!= k)
+    if(i1 != 1 && TKresult[nres]!= k)
       continue;
     if(invalidvarcomb[k]){
       printf("\nCombination (%d) projection ignored because no cases \n",k); 
@@ -8244,6 +8501,11 @@ int decoderesult ( char resultline[], in
   if (strlen(resultsav) >1){
     j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' */
   }
+  if(j == 0){ /* Resultline but no = */
+    TKresult[nres]=0; /* Combination for the nresult and the model */
+    return (0);
+  }
+    
   if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */
     printf("ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs);
     fprintf(ficlog,"ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs);
@@ -8843,7 +9105,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
 }
 
 int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )
-{
+{/* Check ages at death */
   int i, m;
   int firstone=0;
   
@@ -8858,16 +9120,16 @@ int calandcheckages(int imx, int maxwav,
 	*nberr = *nberr + 1;
 	if(firstone == 0){
 	  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);
+	printf("Warning (#%d)! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown but status is a death state %d at wave %d. If you don't know the vital status, please enter -2. If he/she is still alive but don't know the state, please code with '-1 or '.'. Here, we do not believe in a death, skipped.\nOther similar cases in log file\n", *nberr,(int)moisdc[i],(int)andc[i],num[i],i,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;
+	fprintf(ficlog,"Warning (#%d)! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown but status is a death state %d at wave %d. If you don't know the vital status, please enter -2. If he/she is still alive but don't know the state, please code with '-1 or '.'. Here, we do not believe in a death, skipped.\nOther similar cases in log file\n", *nberr,(int)moisdc[i],(int)andc[i],num[i],i,s[m][i],m);
+	s[m][i]=-1;  /* Droping the death status */
       }
       if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
 	(*nberr)++;
-	printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); 
-	fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); 
-	s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */
+	printf("Error (#%d)! Month of death of individual %ld on line %d was unknown (%2d) (year of death is %4d) and status is a death state %d at wave %d. Please impute an arbitrary (or not) month and rerun. Currently this transition to death will be skipped (status is set to -2).\nOther similar cases in log file\n", *nberr, num[i],i,(int)moisdc[i],(int)andc[i],s[m][i],m);
+	fprintf(ficlog,"Error (#%d)! Month of death of individual %ld on line %d was unknown (%2d) (year of death is %4d) and status is a death state %d at wave %d. Please impute an arbitrary (or not) month and rerun. Currently this transition to death will be skipped (status is set to -2).\nOther similar cases in log file\n", *nberr, num[i],i,(int)moisdc[i],(int)andc[i],s[m][i],m);
+	s[m][i]=-2; /* We prefer to skip it (and to skip it in version 0.8a1 too */
       }
     }
   }
@@ -9180,7 +9442,7 @@ int prevalence_limit(double *p, double *
 
   for(k=1; k<=i1;k++){ /* For each combination k of dummy covariates in the model */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
 	continue;
 
       /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
@@ -9277,7 +9539,7 @@ int back_prevalence_limit(double *p, dou
   
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-      if(TKresult[nres]!= k)
+     if(i1 != 1 && TKresult[nres]!= k)
 	continue;
       //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
       fprintf(ficresplb,"#******");
@@ -9336,6 +9598,7 @@ int back_prevalence_limit(double *p, dou
 	fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
       } /* Age */
       /* was end of cptcod */
+      /*fprintf(ficresplb,"\n");*/ /* Seems to be necessary for gnuplot only if two result lines and no covariate. */
     } /* end of any combination */
   } /* end of nres */  
   /* hBijx(p, bage, fage); */
@@ -9380,7 +9643,7 @@ int hPijx(double *p, int bage, int fage)
 		/*  	k=k+1;  */
     for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k=1; k<=i1;k++){
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
 	continue;
       fprintf(ficrespij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++) 
@@ -9452,14 +9715,14 @@ int hPijx(double *p, int bage, int fage)
   
   /* 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 ");
+  fprintf(ficrespijb,"#****** h Bij x Back probability to be in state i at age x-h being in j at x: B1j+B2j+...=1 ");
   i1= pow(2,cptcoveff);
   /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
   /*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
   /*  	k=k+1;  */
   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
 	continue;
       fprintf(ficrespijb,"\n#****** ");
       for(j=1;j<=cptcoveff;j++)
@@ -9486,7 +9749,7 @@ int hPijx(double *p, int bage, int fage)
 	/* 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=");
+	fprintf(ficrespijb,"# Cov Agex agex-h hbijx with i,j=");
 	for(i=1; i<=nlstate;i++)
 	  for(j=1; j<=nlstate+ndeath;j++)
 	    fprintf(ficrespijb," %1d-%1d",i,j);
@@ -9532,6 +9795,7 @@ int main(int argc, char *argv[])
   int NDIM=2;
   int vpopbased=0;
   int nres=0;
+  int endishere=0;
 
   char ca[32], cb[32];
   /*  FILE *fichtm; *//* Html File */
@@ -9576,7 +9840,8 @@ int main(int argc, char *argv[])
   double **prlim;
   double **bprlim;
   double ***param; /* Matrix of parameters */
-  double  *p;
+  double ***paramstart; /* Matrix of starting parameter values */
+  double  *p, *pstart; /* p=param[1][1] pstart is for starting values guessed by freqsummary */
   double **matcov; /* Matrix of covariance */
   double **hess; /* Hessian matrix */
   double ***delti3; /* Scale */
@@ -9922,6 +10187,7 @@ int main(int argc, char *argv[])
     ungetc(c,ficpar);
     
     param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
+    paramstart= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
     for(i=1; i <=nlstate; i++){
       j=0;
       for(jj=1; jj <=nlstate+ndeath; jj++){
@@ -9958,8 +10224,9 @@ run imach with mle=-1 to get a correct t
     }  
     fflush(ficlog);
     
-    /* Reads scales values */
+    /* Reads parameters values */
     p=param[1][1];
+    pstart=paramstart[1][1];
     
     /* Reads comments: lines beginning with '#' */
     while((c=getc(ficpar))=='#' && c!= EOF){
@@ -10224,7 +10491,7 @@ 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);
@@ -10405,7 +10672,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
   /* Calculates basic frequencies. Computes observed prevalence at single age 
 		 and for any valid combination of covariates
      and prints on file fileres'p'. */
-  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \
+  freqsummary(fileres, p, pstart, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \
 	      firstpass, lastpass,  stepm,  weightopt, model);
 
   fprintf(fichtm,"\n");
@@ -10427,9 +10694,9 @@ Interval (in months) between two waves:
   /* For mortality only */
   if (mle==-3){
     ximort=matrix(1,NDIM,1,NDIM); 
-		for(i=1;i<=NDIM;i++)
-			for(j=1;j<=NDIM;j++)
-				ximort[i][j]=0.;
+    for(i=1;i<=NDIM;i++)
+      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);
@@ -10884,70 +11151,31 @@ Please run with mle=-1 to get a correct
     fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d, ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
 		
     /* Other stuffs, more or less useful */    
-    while((c=getc(ficpar))=='#' && c!= EOF){
-      ungetc(c,ficpar);
-      fgets(line, MAXLINE, ficpar);
-      fputs(line,stdout);
-      fputs(line,ficparo);
-    }
-    ungetc(c,ficpar);
-    
-    fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
-    fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
-    fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
-    printf("begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
-    fprintf(ficlog,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
-    
-    while((c=getc(ficpar))=='#' && c!= EOF){
-      ungetc(c,ficpar);
-      fgets(line, MAXLINE, ficpar);
-      fputs(line,stdout);
-      fputs(line,ficparo);
-    }
-    ungetc(c,ficpar);
-    
-    
-    dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.;
-    dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.;
-    
-    fscanf(ficpar,"pop_based=%d\n",&popbased);
-    fprintf(ficlog,"pop_based=%d\n",popbased);
-    fprintf(ficparo,"pop_based=%d\n",popbased);   
-    fprintf(ficres,"pop_based=%d\n",popbased);   
-    
-    while((c=getc(ficpar))=='#' && c!= EOF){
-      ungetc(c,ficpar);
-      fgets(line, MAXLINE, ficpar);
-      fputs(line,stdout);
-      fputs(line,ficres);
-      fputs(line,ficparo);
+    while(fgets(line, MAXLINE, ficpar)) {
+      /* If line starts with a # it is a comment */
+      if (line[0] == '#') {
+	numlinepar++;
+	fputs(line,stdout);
+	fputs(line,ficparo);
+	fputs(line,ficlog);
+	continue;
+      }else
+	break;
     }
-    ungetc(c,ficpar);
-    
-    fscanf(ficpar,"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);
-    fprintf(ficparo,"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);
-    printf("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);
-    fprintf(ficlog,"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);
-    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);
-      fputs(line,ficres);
+
+    if((num_filled=sscanf(line,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav)) !=EOF){
+      
+      if (num_filled != 7) {
+	printf("Error: Not 7 (data)parameters in line but %d, for example:begin-prev-date=1/1/1990 end-prev-date=1/6/2004  mov_average=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	fprintf(ficlog,"Error: Not 7 (data)parameters in line but %d, for example:begin-prev-date=1/1/1990 end-prev-date=1/6/2004  mov_average=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	goto end;
+      }
+      printf("begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
+      fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
+      fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
+      fprintf(ficlog,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
     }
-    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.*/
-    
-    /* Results */
-    nresult=0;
+
     while(fgets(line, MAXLINE, ficpar)) {
       /* If line starts with a # it is a comment */
       if (line[0] == '#') {
@@ -10955,49 +11183,114 @@ Please run with mle=-1 to get a correct
 	fputs(line,stdout);
 	fputs(line,ficparo);
 	fputs(line,ficlog);
-	fputs(line,ficres);
 	continue;
       }else
 	break;
     }
-    if (!feof(ficpar))
-    while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){
-      if (num_filled == 0){
-	resultline[0]='\0';
-      break;
-      } else if (num_filled != 1){
-	printf("ERROR %d: result line should be at minimum 'result=' %s\n",num_filled, line);
-      }
-      nresult++; /* Sum of resultlines */
-      printf("Result %d: result=%s\n",nresult, resultline);
-      if(nresult > MAXRESULTLINES){
-	printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult);
-	fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult);
+    
+    
+    dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.;
+    dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.;
+    
+    if((num_filled=sscanf(line,"pop_based=%d\n",&popbased)) !=EOF){
+      if (num_filled != 1) {
+	printf("Error: Not 1 (data)parameters in line but %d, for example:pop_based=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	fprintf(ficlog,"Error: Not 1 (data)parameters in line but %d, for example: pop_based=1\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
 	goto end;
       }
-      decoderesult(resultline, nresult); /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */
-      fprintf(ficparo,"result: %s\n",resultline);
-      fprintf(ficres,"result: %s\n",resultline);
-      fprintf(ficlog,"result: %s\n",resultline);
-      while(fgets(line, MAXLINE, ficpar)) {
+      printf("pop_based=%d\n",popbased);
+      fprintf(ficlog,"pop_based=%d\n",popbased);
+      fprintf(ficparo,"pop_based=%d\n",popbased);   
+      fprintf(ficres,"pop_based=%d\n",popbased);   
+    }
+     
+    /* Results */
+    nresult=0;
+    do{
+      if(!fgets(line, MAXLINE, ficpar)){
+	endishere=1;
+	parameterline=14;
+      }else if (line[0] == '#') {
 	/* If line starts with a # it is a comment */
-	if (line[0] == '#') {
-	  numlinepar++;
-	  fputs(line,stdout);
-	  fputs(line,ficparo);
-	  fputs(line,ficres);
-	  fputs(line,ficlog);
-	  continue;
-	}else
-	  break;
+	numlinepar++;
+	fputs(line,stdout);
+	fputs(line,ficparo);
+	fputs(line,ficlog);
+	continue;
+      }else if(sscanf(line,"prevforecast=%[^\n]\n",modeltemp))
+	parameterline=11;
+      else if(sscanf(line,"backcast=%[^\n]\n",modeltemp))
+	parameterline=12;
+      else if(sscanf(line,"result:%[^\n]\n",modeltemp))
+	parameterline=13;
+      else{
+	parameterline=14;
       }
-      if (feof(ficpar))
+      switch (parameterline){ 
+      case 11:
+	if((num_filled=sscanf(line,"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)) !=EOF){
+	  if (num_filled != 8) {
+	    printf("Error: Not 8 (data)parameters in line but %d, for example:prevforecast=1 starting-proj-date=1/1/1990 final-proj-date=1/1/2000 mobil_average=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	    fprintf(ficlog,"Error: Not 8 (data)parameters in line but %d, for example:prevforecast=1 starting-proj-date=1/1/1990 final-proj-date=1/1/2000 mov_average=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	    goto end;
+	  }
+	  fprintf(ficparo,"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);
+	  printf("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);
+	  fprintf(ficlog,"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);
+	  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.*/
+	}
 	break;
-      else{ /* Processess output results for this combination of covariate values */
-      }				   
-    } /* end while */
-
-
+      case 12:
+	/*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);*/
+	if((num_filled=sscanf(line,"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)) !=EOF){
+	  if (num_filled != 8) {
+	    printf("Error: Not 8 (data)parameters in line but %d, for example:backcast=1 starting-back-date=1/1/1990 finloal-back-date=1/1/1970 mobil_average=1\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	    fprintf(ficlog,"Error: Not 8 (data)parameters in line but %d, for example:backcast=1 starting-back-date=1/1/1990 finloal-back-date=1/1/1970 mobil_average=1\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
+	    goto end;
+	  }
+	  printf("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.*/
+	}
+	break;
+      case 13:
+	if((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){
+	  if (num_filled == 0){
+	    resultline[0]='\0';
+	    printf("Warning %d: no result line! It should be at minimum 'result: V2=0 V1=1 or result:.\n%s\n", num_filled, line);
+	    fprintf(ficlog,"Warning %d: no result line! It should be at minimum 'result: V2=0 V1=1 or result:.\n%s\n", num_filled, line);
+	    break;
+	  } else if (num_filled != 1){
+	    printf("ERROR %d: result line! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",num_filled, line);
+	    fprintf(ficlog,"ERROR %d: result line! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",num_filled, line);
+	  }
+	  nresult++; /* Sum of resultlines */
+	  printf("Result %d: result=%s\n",nresult, resultline);
+	  if(nresult > MAXRESULTLINES){
+	    printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult);
+	    fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult);
+	    goto end;
+	  }
+	  decoderesult(resultline, nresult); /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */
+	  fprintf(ficparo,"result: %s\n",resultline);
+	  fprintf(ficres,"result: %s\n",resultline);
+	  fprintf(ficlog,"result: %s\n",resultline);
+	  break;
+	case 14: 
+	  if(ncovmodel >2 && nresult==0 ){
+	    printf("ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line);
+	    goto end;
+	  }
+	  break;
+	default:
+	  nresult=1;
+	  decoderesult(".",nresult ); /* No covariate */
+	}
+      } /* End switch parameterline */
+    }while(endishere==0); /* End do */
     
     /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */
     /* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */
@@ -11014,7 +11307,7 @@ Please run with mle=-1 to get a correct
       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,backcast, estepm, \
+		 model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \
 		 jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2);
 		
     /*------------ free_vector  -------------*/
@@ -11067,6 +11360,7 @@ Please run with mle=-1 to get a correct
       mobaverage=mobaverages;
       if (mobilav!=0) {
 	printf("Movingaveraging observed prevalence\n");
+	fprintf(ficlog,"Movingaveraging observed prevalence\n");
 	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);
@@ -11076,6 +11370,7 @@ Please run with mle=-1 to get a correct
       /* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
       else if (mobilavproj !=0) {
 	printf("Movingaveraging projected observed prevalence\n");
+	fprintf(ficlog,"Movingaveraging projected observed prevalence\n");
 	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);
@@ -11138,7 +11433,7 @@ Please run with mle=-1 to get a correct
     
     for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
 	continue;
       fprintf(ficreseij,"\n#****** ");
       printf("\n#****** ");
@@ -11211,7 +11506,7 @@ Please run with mle=-1 to get a correct
     
     for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
 	continue;
       printf("\n#****** Result for:");
       fprintf(ficrest,"\n#****** Result for:");
@@ -11352,7 +11647,7 @@ Please run with mle=-1 to get a correct
 
     for(nres=1; nres <= nresult; nres++) /* For each resultline */
     for(k=1; k<=i1;k++){
-      if(TKresult[nres]!= k)
+      if(i1 != 1 && TKresult[nres]!= k)
 	continue;
       fprintf(ficresvpl,"\n#****** ");
       printf("\n#****** ");