--- imach/src/imach.c	2017/06/30 13:39:33	1.275
+++ imach/src/imach.c	2018/04/21 21:02:16	1.285
@@ -1,6 +1,38 @@
-/* $Id: imach.c,v 1.275 2017/06/30 13:39:33 brouard Exp $
+/* $Id: imach.c,v 1.285 2018/04/21 21:02:16 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.285  2018/04/21 21:02:16  brouard
+  Summary: Some bugs fixed, valgrind tested
+
+  Revision 1.284  2018/04/20 05:22:13  brouard
+  Summary: Computing mean and stdeviation of fixed quantitative variables
+
+  Revision 1.283  2018/04/19 14:49:16  brouard
+  Summary: Some minor bugs fixed
+
+  Revision 1.282  2018/02/27 22:50:02  brouard
+  *** empty log message ***
+
+  Revision 1.281  2018/02/27 19:25:23  brouard
+  Summary: Adding second argument for quitting
+
+  Revision 1.280  2018/02/21 07:58:13  brouard
+  Summary: 0.99r15
+
+  New Makefile with recent VirtualBox 5.26. Bug in sqrt negatve in imach.c
+
+  Revision 1.279  2017/07/20 13:35:01  brouard
+  Summary: temporary working
+
+  Revision 1.278  2017/07/19 14:09:02  brouard
+  Summary: Bug for mobil_average=0 and prevforecast fixed(?)
+
+  Revision 1.277  2017/07/17 08:53:49  brouard
+  Summary: BOM files can be read now
+
+  Revision 1.276  2017/06/30 15:48:31  brouard
+  Summary: Graphs improvements
+
   Revision 1.275  2017/06/30 13:39:33  brouard
   Summary: Saito's color
 
@@ -1028,12 +1060,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.275 2017/06/30 13:39:33 brouard Exp $ */
+/* $Id: imach.c,v 1.285 2018/04/21 21:02:16 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.275 $ $Date: 2017/06/30 13:39:33 $"; 
+char copyright[]="April 2018,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018";
+char fullversion[]="$Revision: 1.285 $ $Date: 2018/04/21 21:02:16 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -2507,15 +2539,18 @@ void powell(double p[], double **xi, int
   
   double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij, int nres)
   {
-    /* Computes the prevalence limit in each live state at age x and for covariate combination ij 
-       (and selected quantitative values in nres)
-       by left multiplying the unit
-       matrix by transitions matrix until convergence is reached with precision ftolpl */
-  /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1  = Wx-n Px-n ... Px-2 Px-1 I */
-  /* Wx is row vector: population in state 1, population in state 2, population dead */
-  /* or prevalence in state 1, prevalence in state 2, 0 */
-  /* newm is the matrix after multiplications, its rows are identical at a factor */
-  /* Initial matrix pimij */
+    /**< Computes the prevalence limit in each live state at age x and for covariate combination ij 
+     *   (and selected quantitative values in nres)
+     *  by left multiplying the unit
+     *  matrix by transitions matrix until convergence is reached with precision ftolpl 
+     * Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1  = Wx-n Px-n ... Px-2 Px-1 I
+     * Wx is row vector: population in state 1, population in state 2, population dead
+     * or prevalence in state 1, prevalence in state 2, 0
+     * newm is the matrix after multiplications, its rows are identical at a factor.
+     * Inputs are the parameter, age, a tolerance for the prevalence limit ftolpl.
+     * Output is prlim.
+     * Initial matrix pimij 
+     */
   /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */
   /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */
   /*  0,                   0                  , 1} */
@@ -4342,7 +4377,7 @@ void  freqsummary(char fileres[], double
   double ***freq; /* Frequencies */
   double *x, *y, a=0.,b=0.,r=1., sa=0., sb=0.; /* for regression, y=b+m*x and r is the correlation coefficient */
   int no=0, linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb);
-  double *meanq;
+  double *meanq, *stdq, *idq;
   double **meanqt;
   double *pp, **prop, *posprop, *pospropt;
   double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0;
@@ -4355,6 +4390,8 @@ void  freqsummary(char fileres[], double
   pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ 
   /* prop=matrix(1,nlstate,iagemin,iagemax+3); */
   meanq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */
+  stdq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */
+  idq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */
   meanqt=matrix(1,lastpass,1,nqtveff);
   strcpy(fileresp,"P_");
   strcat(fileresp,fileresu);
@@ -4463,13 +4500,16 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	posprop[i]=0;
 	pospropt[i]=0;
       }
-      /* for (z1=1; z1<= nqfveff; z1++) {   */
-      /*   meanq[z1]+=0.; */
+      for (z1=1; z1<= nqfveff; z1++) { /* zeroing for each combination j1 as well as for the total */
+        idq[z1]=0.;
+        meanq[z1]=0.;
+        stdq[z1]=0.;
+      }
+      /* for (z1=1; z1<= nqtveff; z1++) { */
       /*   for(m=1;m<=lastpass;m++){ */
-      /* 	meanqt[m][z1]=0.; */
-      /*   } */
-      /* } */
-      
+      /* 	  meanqt[m][z1]=0.; */
+      /* 	} */
+      /* }       */
       /* dateintsum=0; */
       /* k2cpt=0; */
       
@@ -4479,9 +4519,6 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	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]; *\/ */
@@ -4502,7 +4539,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	}/* end j==0 */
 	if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */
 	  /* for(m=firstpass; m<=lastpass; m++){ */
-	  for(mi=1; mi<wav[iind];mi++){ /* For that wave */
+	  for(mi=1; mi<wav[iind];mi++){ /* For each wave */
 	    m=mw[mi][iind];
 	    if(j!=0){
 	      if(anyvaryingduminmodel==1){ /* Some are varying covariates */
@@ -4522,7 +4559,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	      }/* 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){
+	    if(bool==1){ /*Selected */
 	      /* 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*/
@@ -4540,6 +4577,11 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 		  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 */
+		  for (z1=1; z1<= nqfveff; z1++) { /* Quantitative variables, calculating mean */
+		    idq[z1]=idq[z1]+weight[iind];
+		    meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind];  /* Computes mean of quantitative with selected filter */
+		    stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /* *weight[iind];*/  /* Computes mean of quantitative with selected filter */
+		  }
 		  /* 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]; */
@@ -4555,6 +4597,11 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	      bool=1;
 	    }/* end bool 2 */
 	  } /* end m */
+	  /* for (z1=1; z1<= nqfveff; z1++) { /\* Quantitative variables, calculating mean *\/ */
+	  /*   idq[z1]=idq[z1]+weight[iind]; */
+	  /*   meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind];  /\* Computes mean of quantitative with selected filter *\/ */
+	  /*   stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /\* *weight[iind];*\/  /\* Computes mean of quantitative with selected filter *\/ */
+	  /* } */
 	} /* end bool */
       } /* end iind = 1 to imx */
       /* prop[s][age] is feeded for any initial and valid live state as well as
@@ -4592,6 +4639,27 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	fprintf(ficresphtmfr, "**********</h3>\n");
 	fprintf(ficlog, "**********\n");
       }
+      /*
+	Printing means of quantitative variables if any
+      */
+      for (z1=1; z1<= nqfveff; z1++) {
+	fprintf(ficlog,"Mean of fixed quantitative variable V%d on %.0f individuals sum=%f", ncovcol+z1, idq[z1], meanq[z1]);
+	fprintf(ficlog,", mean=%.3g\n",meanq[z1]/idq[z1]);
+	if(weightopt==1){
+	  printf(" Weighted mean and standard deviation of");
+	  fprintf(ficlog," Weighted mean and standard deviation of");
+	  fprintf(ficresphtmfr," Weighted mean and standard deviation of");
+	}
+	printf(" fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1]));
+	fprintf(ficlog," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1]));
+	fprintf(ficresphtmfr," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)<p>\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1]));
+      }
+      /* for (z1=1; z1<= nqtveff; z1++) { */
+      /* 	for(m=1;m<=lastpass;m++){ */
+      /* 	  fprintf(ficresphtmfr,"V quantitative id %d, pass id=%d, mean=%f<p>\n", z1, m, meanqt[m][z1]); */
+      /*   } */
+      /* } */
+
       fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
       if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */
         fprintf(ficresp, " Age");
@@ -4826,7 +4894,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
 	    fprintf(ficlog,"\n");
 	  }
 	}
-      }
+      } /* end of state i */
       printf("#Freqsummary\n");
       fprintf(ficlog,"\n");
       for(s1=-1; s1 <=nlstate+ndeath; s1++){
@@ -4872,7 +4940,9 @@ Title=%s <br>Datafile=%s Firstpass=%d La
   fclose(ficresp);
   fclose(ficresphtm);
   fclose(ficresphtmfr);
+  free_vector(idq,1,nqfveff);
   free_vector(meanq,1,nqfveff);
+  free_vector(stdq,1,nqfveff);
   free_matrix(meanqt,1,lastpass,1,nqtveff);
   free_vector(x, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
   free_vector(y, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
@@ -5288,6 +5358,9 @@ void  concatwav(int wav[], int **dh, int
    /* *cptcov=0; */
  
    for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */
+   for (k=1; k <= maxncov; k++)
+     for(j=1; j<=2; j++)
+       nbcode[k][j]=0; /* Valgrind */
 
    /* Loop on covariates without age and products and no quantitative variable */
    /* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */
@@ -5763,10 +5836,11 @@ void  concatwav(int wav[], int **dh, int
 /************ Variance ******************/
  void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[], int nres)
  {
-   /* Variance of health expectancies */
-   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
-   /* double **newm;*/
-   /* int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav)*/
+   /** Variance of health expectancies 
+    *  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);
+    * double **newm;
+    * int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav) 
+    */
   
    /* int movingaverage(); */
    double **dnewm,**doldm;
@@ -5774,11 +5848,11 @@ void  concatwav(int wav[], int **dh, int
    int i, j, nhstepm, hstepm, h, nstepm ;
    int k;
    double *xp;
-   double **gp, **gm;  /* for var eij */
-   double ***gradg, ***trgradg; /*for var eij */
-   double **gradgp, **trgradgp; /* for var p point j */
-   double *gpp, *gmp; /* for var p point j */
-   double **varppt; /* for var p point j nlstate to nlstate+ndeath */
+   double **gp, **gm;  /**< for var eij */
+   double ***gradg, ***trgradg; /**< for var eij */
+   double **gradgp, **trgradgp; /**< for var p point j */
+   double *gpp, *gmp; /**< for var p point j */
+   double **varppt; /**< for var p point j nlstate to nlstate+ndeath */
    double ***p3mat;
    double age,agelim, hf;
    /* double ***mobaverage; */
@@ -5839,7 +5913,7 @@ void  concatwav(int wav[], int **dh, int
    /* fprintf(fichtm, "#Local time at start: %s", strstart);*/
    fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
    fprintf(fichtm,"\n<br>%s  <br>\n",digitp);
-   /*   } */
+
    varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
    pstamp(ficresvij);
    fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are ");
@@ -5894,9 +5968,12 @@ void  concatwav(int wav[], int **dh, int
        for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/
 	 xp[i] = x[i] + (i==theta ?delti[theta]:0);
        }
-			
+       /**< Computes the prevalence limit with parameter theta shifted of delta up to ftolpl precision and 
+	* returns into prlim .
+	*/		
        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres);
-			
+
+       /* If popbased = 1 we use crossection prevalences. Previous step is useless but prlim is created */
        if (popbased==1) {
 	 if(mobilav ==0){
 	   for(i=1; i<=nlstate;i++)
@@ -5906,23 +5983,28 @@ void  concatwav(int wav[], int **dh, int
 	     prlim[i][i]=mobaverage[(int)age][i][ij];
 	 }
        }
-			
-       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres);  /* Returns p3mat[i][j][h] for h=1 to nhstepm */
+       /**< Computes the shifted transition matrix \f$ {}{h}_p^{ij}_x\f$ at horizon h.
+	*/			
+       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres);  /* Returns p3mat[i][j][h] for h=0 to nhstepm */
+       /**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}_x\f$, which are the probability
+	* at horizon h in state j including mortality.
+	*/
        for(j=1; j<= nlstate; j++){
 	 for(h=0; h<=nhstepm; h++){
 	   for(i=1, gp[h][j]=0.;i<=nlstate;i++)
 	     gp[h][j] += prlim[i][i]*p3mat[i][j][h];
 	 }
        }
-       /* Next for computing probability of death (h=1 means
+       /* Next for computing shifted+ probability of death (h=1 means
 	  computed over hstepm matrices product = hstepm*stepm months) 
-	  as a weighted average of prlim.
+	  as a weighted average of prlim(i) * p(i,j) p.3=w1*p13 + w2*p23 .
        */
        for(j=nlstate+1;j<=nlstate+ndeath;j++){
 	 for(i=1,gpp[j]=0.; i<= nlstate; i++)
 	   gpp[j] += prlim[i][i]*p3mat[i][j][1];
-       }    
-       /* end probability of death */
+       }
+       
+       /* Again with minus shift */
 			
        for(i=1; i<=npar; i++) /* Computes gradient x - delta */
 	 xp[i] = x[i] - (i==theta ?delti[theta]:0);
@@ -5955,19 +6037,23 @@ void  concatwav(int wav[], int **dh, int
 	 for(i=1,gmp[j]=0.; i<= nlstate; i++)
 	   gmp[j] += prlim[i][i]*p3mat[i][j][1];
        }    
-       /* end probability of death */
-			
+       /* end shifting computations */
+
+       /**< Computing gradient matrix at horizon h 
+	*/
        for(j=1; j<= nlstate; j++) /* vareij */
 	 for(h=0; h<=nhstepm; h++){
 	   gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
 	 }
-			
-       for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */
+       /**< Gradient of overall mortality p.3 (or p.j) 
+	*/
+       for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu mortality from j */
 	 gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
        }
 			
      } /* End theta */
-		
+     
+     /* We got the gradient matrix for each theta and state j */		
      trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */
 		
      for(h=0; h<=nhstepm; h++) /* veij */
@@ -5978,13 +6064,19 @@ void  concatwav(int wav[], int **dh, int
      for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */
        for(theta=1; theta <=npar; theta++)
 	 trgradgp[j][theta]=gradgp[theta][j];
-		
+     /**< as well as its transposed matrix 
+      */		
 		
      hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
      for(i=1;i<=nlstate;i++)
        for(j=1;j<=nlstate;j++)
 	 vareij[i][j][(int)age] =0.;
-		
+
+     /* Computing trgradg by matcov by gradg at age and summing over h
+      * and k (nhstepm) formula 15 of article
+      * Lievre-Brouard-Heathcote
+      */
+     
      for(h=0;h<=nhstepm;h++){
        for(k=0;k<=nhstepm;k++){
 	 matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);
@@ -5995,7 +6087,11 @@ void  concatwav(int wav[], int **dh, int
        }
      }
 		
-     /* pptj */
+     /* pptj is p.3 or p.j = trgradgp by cov by gradgp, variance of
+      * p.j overall mortality formula 49 but computed directly because
+      * we compute the grad (wix pijx) instead of grad (pijx),even if
+      * wix is independent of theta.
+      */
      matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);
      matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);
      for(j=nlstate+1;j<=nlstate+ndeath;j++)
@@ -6604,7 +6700,12 @@ To be simple, these graphs help to under
 		 }
 								
 		 /* Eigen vectors */
-		 v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
+		 if(1+(v1-lc1)*(v1-lc1)/cv12/cv12 <1.e-5){
+		   printf(" Error sqrt of a negative number: %lf\n",1+(v1-lc1)*(v1-lc1)/cv12/cv12);
+		   fprintf(ficlog," Error sqrt of a negative number: %lf\n",1+(v1-lc1)*(v1-lc1)/cv12/cv12);
+		   v11=(1./sqrt(fabs(1+(v1-lc1)*(v1-lc1)/cv12/cv12)));
+		 }else
+		   v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
 		 /*v21=sqrt(1.-v11*v11); *//* error */
 		 v21=(lc1-v1)/cv12*v11;
 		 v12=-v21;
@@ -6635,8 +6736,8 @@ To be simple, these graphs help to under
 		   fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
 		   fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
 		   fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",	\
-			   mu1,std,v11,sqrt(lc1),v12,sqrt(fabs(lc2)),	\
-			   mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(lc2))); /* For gnuplot only */
+			   mu1,std,v11,sqrt(fabs(lc1)),v12,sqrt(fabs(lc2)), \
+			   mu2,std,v21,sqrt(fabs(lc1)),v22,sqrt(fabs(lc2))); /* For gnuplot only */
 		 }else{
 		   first=0;
 		   fprintf(fichtmcov," %d (%.3f),",(int) age, c12);
@@ -6811,7 +6912,7 @@ divided by h: <sub>h</sub>P<sub>ij</sub>
      for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\
  Or probability to survive in various states (1 to %d) being in state %d at different ages.	\
- <a href=\"%s_%d-%d-%d.svg\">%s_%d%d-%d.svg</a><br> <img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
+ <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> <img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
      }
      /* Period (stable) prevalence in each health state */
      for(cpt=1; cpt<=nlstate;cpt++){
@@ -7052,7 +7153,8 @@ void printinggnuplot(char fileresu[], ch
       
 	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 label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
+	/* fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); */
+	fprintf(ficgp,"set title \"Alive state %d %s\" font \"Helvetica,12\"\n",cpt,gplotlabel);
 	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*/
@@ -7137,7 +7239,7 @@ void printinggnuplot(char fileresu[], ch
 	      if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
 	      else fprintf(ficgp," %%*lf (%%*lf)");
 	    } 
-	    fprintf(ficgp,"\" t\"95%% CI\" w l lt 5,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres); 
+	    fprintf(ficgp,"\" t\"95%% CI\" w l lt 4,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres); 
 	    for (i=1; i<= nlstate ; i ++) {
 	      if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
 	      else fprintf(ficgp," %%*lf (%%*lf)");
@@ -7145,7 +7247,8 @@ void printinggnuplot(char fileresu[], ch
 	    fprintf(ficgp,"\" t\"\" w l lt 4");
 	  } /* end if backprojcast */
 	} /* end if backcast */
-	fprintf(ficgp,"\nset out ;unset label;\n");
+	/* fprintf(ficgp,"\nset out ;unset label;\n"); */
+	fprintf(ficgp,"\nset out ;unset title;\n");
       } /* nres */
     } /* k1 */
   } /* cpt */
@@ -7759,7 +7862,7 @@ set ter svg size 640, 480\nunset log y\n
 	continue;
       fprintf(ficgp,"\n\n# Combination of dummy  k1=%d which is ",k1);
       strcpy(gplotlabel,"(");
-      sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1);
+      /*sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1);*/
       for (k=1; k<=cptcoveff; k++){    /* For each correspondig covariate value  */
 	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
 	/* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
@@ -7776,7 +7879,9 @@ set ter svg size 640, 480\nunset log y\n
       strcpy(gplotlabel+strlen(gplotlabel),")");
       fprintf(ficgp,"\n#\n");
       fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),k1,ng,nres);
-      fprintf(ficgp,"\nset label \"%s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",gplotlabel);
+      fprintf(ficgp,"\nset key outside ");
+      /* fprintf(ficgp,"\nset label \"%s\" at graph 1.2,0.5 center rotate font \"Helvetica,12\"\n",gplotlabel); */
+      fprintf(ficgp,"\nset title \"%s\" font \"Helvetica,12\"\n",gplotlabel);
       fprintf(ficgp,"\nset ter svg size 640, 480 ");
       if (ng==1){
 	fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */
@@ -7896,12 +8001,12 @@ set ter svg size 640, 480\nunset log y\n
 	    }
 	    fprintf(ficgp,")");
 	    if(ng ==2)
-	      fprintf(ficgp," w l lt (%d*%d+%d)%%%d+1 dt %d t \"p%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
+	      fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"p%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
 	    else /* ng= 3 */
-	      fprintf(ficgp," w l lt (%d*%d+%d)%%%d+1 dt %d t \"i%d%d\" ",  nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
+	      fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"i%d%d\" ",  nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
 	  }else{ /* end ng <> 1 */
 	    if( k !=k2) /* logit p11 is hard to draw */
-	      fprintf(ficgp," w l lt (%d*%d+%d)%%%d+1 dt %d t \"logit(p%d%d)\" ",  nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
+	      fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"logit(p%d%d)\" ",  nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
 	  }
 	  if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)
 	    fprintf(ficgp,",");
@@ -7910,7 +8015,8 @@ set ter svg size 640, 480\nunset log y\n
 	  i=i+ncovmodel;
 	} /* end k */
       } /* end k2 */
-      fprintf(ficgp,"\n set out; unset label;\n");
+      /* fprintf(ficgp,"\n set out; unset label;set key default;\n"); */
+      fprintf(ficgp,"\n set out; unset title;set key default;\n");
     } /* end k1 */
   } /* end ng */
   /* avoid: */
@@ -7934,8 +8040,8 @@ set ter svg size 640, 480\nunset log y\n
    double *agemingoodr, *agemaxgoodr; 
   
   
-   /* modcovmax=2*cptcoveff;/\* Max number of modalities. We suppose  */
-   /* 		   a covariate has 2 modalities, should be equal to ncovcombmax  *\/ */
+   /* modcovmax=2*cptcoveff;  Max number of modalities. We suppose  */
+   /* 		   a covariate has 2 modalities, should be equal to ncovcombmax   */
 
    sumnewp = vector(1,ncovcombmax);
    sumnewm = vector(1,ncovcombmax);
@@ -8259,11 +8365,11 @@ set ter svg size 640, 480\nunset log y\n
 	for(j=1; j<=nlstate+ndeath;j++) {
 	  ppij=0.;
 	  for(i=1; i<=nlstate;i++) {
-	    /* if (mobilav>=1)  */
-	    ppij=ppij+p3mat[i][j][h]*prev[(int)agec][i][k];
-	    /* else { */ /* even if mobilav==-1 we use mobaverage */
-	    /* 	ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
-	    /* } */
+	    if (mobilav>=1)
+	     ppij=ppij+p3mat[i][j][h]*prev[(int)agec][i][k];
+	    else { /* even if mobilav==-1 we use mobaverage, probs may not sums to 1 */
+	    	ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k];
+	    }
 	    fprintf(ficresf," %.3f", p3mat[i][j][h]);
 	  } /* end i */
 	  fprintf(ficresf," %.3f", ppij);
@@ -9630,7 +9736,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
 Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for  product \n\
 Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
 Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
-  for(k=1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;}
+  for(k=-1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;}
   for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */
     if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
       Fixed[k]= 0;
@@ -9880,11 +9986,12 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
   /* Searching for doublons in the model */
   for(k1=1; k1<= cptcovt;k1++){
     for(k2=1; k2 <k1;k2++){
-      if((Typevar[k1]==Typevar[k2]) && (Fixed[Tvar[k1]]==Fixed[Tvar[k2]]) && (Dummy[Tvar[k1]]==Dummy[Tvar[k2]] )){
+      /* if((Typevar[k1]==Typevar[k2]) && (Fixed[Tvar[k1]]==Fixed[Tvar[k2]]) && (Dummy[Tvar[k1]]==Dummy[Tvar[k2]] )){ */
+      if((Typevar[k1]==Typevar[k2]) && (Fixed[k1]==Fixed[k2]) && (Dummy[k1]==Dummy[k2] )){
 	if((Typevar[k1] == 0 || Typevar[k1] == 1)){ /* Simple or age product */
 	  if(Tvar[k1]==Tvar[k2]){
-	    printf("Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
-	    fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
+	    printf("Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[k1],Dummy[k1]);
+	    fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[k1],Dummy[k1]); fflush(ficlog);
 	    return(1);
 	  }
 	}else if (Typevar[k1] ==2){
@@ -10605,7 +10712,7 @@ int main(int argc, char *argv[])
   int vpopbased=0;
   int nres=0;
   int endishere=0;
-
+  int noffset=0;
   int ncurrv=0; /* Temporary variable */
   
   char ca[32], cb[32];
@@ -10742,8 +10849,13 @@ int main(int argc, char *argv[])
       if(pathr[0] == '\0') break; /* Dirty */
     }
   }
+  else if (argc<=2){
+    strcpy(pathtot,argv[1]);
+  }
   else{
     strcpy(pathtot,argv[1]);
+    strcpy(z,argv[2]);
+    printf("\nargv[2]=%s z=%c\n",argv[2],z[0]);
   }
   /*if(getcwd(pathcd, MAXLINE)!= NULL)printf ("Error pathcd\n");*/
   /*cygwin_split_path(pathtot,path,optionfile);
@@ -10821,8 +10933,6 @@ int main(int argc, char *argv[])
     exit(70); 
   }
 
-
-
   strcpy(filereso,"o");
   strcat(filereso,fileresu);
   if((ficparo=fopen(filereso,"w"))==NULL) { /* opened on subdirectory */
@@ -10831,17 +10941,52 @@ int main(int argc, char *argv[])
     fflush(ficlog);
     goto end;
   }
+      /*-------- Rewriting parameter file ----------*/
+  strcpy(rfileres,"r");    /* "Rparameterfile */
+  strcat(rfileres,optionfilefiname);    /* Parameter file first name */
+  strcat(rfileres,".");    /* */
+  strcat(rfileres,optionfilext);    /* Other files have txt extension */
+  if((ficres =fopen(rfileres,"w"))==NULL) {
+    printf("Problem writing new parameter file: %s\n", rfileres);goto end;
+    fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end;
+    fflush(ficlog);
+    goto end;
+  }
+  fprintf(ficres,"#IMaCh %s\n",version);
 
+				      
   /* Reads comments: lines beginning with '#' */
   numlinepar=0;
-
-    /* First parameter line */
+  /* Is it a BOM UTF-8 Windows file? */
+  /* First parameter line */
   while(fgets(line, MAXLINE, ficpar)) {
+    noffset=0;
+    if( line[0] == (char)0xEF && line[1] == (char)0xBB) /* EF BB BF */
+    {
+      noffset=noffset+3;
+      printf("# File is an UTF8 Bom.\n"); // 0xBF
+    }
+    else if( line[0] == (char)0xFE && line[1] == (char)0xFF)
+    {
+      noffset=noffset+2;
+      printf("# File is an UTF16BE BOM file\n");
+    }
+    else if( line[0] == 0 && line[1] == 0)
+    {
+      if( line[2] == (char)0xFE && line[3] == (char)0xFF){
+	noffset=noffset+4;
+	printf("# File is an UTF16BE BOM file\n");
+      }
+    } else{
+      ;/*printf(" Not a BOM file\n");*/
+    }
+  
     /* If line starts with a # it is a comment */
-    if (line[0] == '#') {
+    if (line[noffset] == '#') {
       numlinepar++;
       fputs(line,stdout);
       fputs(line,ficparo);
+      fputs(line,ficres);
       fputs(line,ficlog);
       continue;
     }else
@@ -10851,18 +10996,24 @@ int main(int argc, char *argv[])
 			title, datafile, &lastobs, &firstpass,&lastpass)) !=EOF){
     if (num_filled != 5) {
       printf("Should be 5 parameters\n");
+      fprintf(ficlog,"Should be 5 parameters\n");
     }
     numlinepar++;
     printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass);
+    fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass);
+    fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass);
+    fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass);
   }
   /* Second parameter line */
   while(fgets(line, MAXLINE, ficpar)) {
-    /* If line starts with a # it is a comment */
+    /* while(fscanf(ficpar,"%[^\n]", line)) { */
+    /* If line starts with a # it is a comment. Strangely fgets reads the EOL and fputs doesn't */
     if (line[0] == '#') {
       numlinepar++;
-      fputs(line,stdout);
-      fputs(line,ficparo);
-      fputs(line,ficlog);
+      printf("%s",line);
+      fprintf(ficres,"%s",line);
+      fprintf(ficparo,"%s",line);
+      fprintf(ficlog,"%s",line);
       continue;
     }else
       break;
@@ -10872,8 +11023,13 @@ int main(int argc, char *argv[])
     if (num_filled != 11) {
       printf("Not 11 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nqv=1 ntv=2 nqtv=1  nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n");
       printf("but line=%s\n",line);
+      fprintf(ficlog,"Not 11 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nqv=1 ntv=2 nqtv=1  nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n");
+      fprintf(ficlog,"but line=%s\n",line);
     }
     printf("ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt);
+    fprintf(ficparo,"ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt);
+    fprintf(ficres,"ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt);
+    fprintf(ficlog,"ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt);
   }
   /* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */
   /*ftolpl=6.e-4; *//* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
@@ -10882,22 +11038,18 @@ int main(int argc, char *argv[])
     /* If line starts with a # it is a comment */
     if (line[0] == '#') {
       numlinepar++;
-      fputs(line,stdout);
-      fputs(line,ficparo);
-      fputs(line,ficlog);
+      printf("%s",line);
+      fprintf(ficres,"%s",line);
+      fprintf(ficparo,"%s",line);
+      fprintf(ficlog,"%s",line);
       continue;
     }else
       break;
   }
   if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){
-    if (num_filled == 0){
-      printf("ERROR %d: Model should be at minimum 'model=1+age.' WITHOUT space:'%s'\n",num_filled, line);
-      fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' WITHOUT space:'%s'\n",num_filled, line);
-      model[0]='\0';
-      goto end;
-    } else if (num_filled != 1){
-      printf("ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line);
-      fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line);
+    if (num_filled != 1){
+      printf("ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line);
+      fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line);
       model[0]='\0';
       goto end;
     }
@@ -10910,20 +11062,23 @@ int main(int argc, char *argv[])
     }
     /* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */
     printf("model=1+age+%s\n",model);fflush(stdout);
+    fprintf(ficparo,"model=1+age+%s\n",model);fflush(stdout);
+    fprintf(ficres,"model=1+age+%s\n",model);fflush(stdout);
+    fprintf(ficlog,"model=1+age+%s\n",model);fflush(stdout);
   }
   /* fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */
   /* numlinepar=numlinepar+3; /\* In general *\/ */
   /* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */
-  fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model);
-  fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model);
+  /* fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model); */
+  /* fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model); */
   fflush(ficlog);
   /* if(model[0]=='#'|| model[0]== '\0'){ */
   if(model[0]=='#'){
-    printf("Error in 'model' line: model should start with 'model=1+age+' and end with '.' \n \
- 'model=1+age+.' or 'model=1+age+V1.' or 'model=1+age+age*age+V1+V1*age.' or \n \
- 'model=1+age+V1+V2.' or 'model=1+age+V1+V2+V1*V2.' etc. \n");		\
+    printf("Error in 'model' line: model should start with 'model=1+age+' and end without space \n \
+ 'model=1+age+' or 'model=1+age+V1.' or 'model=1+age+age*age+V1+V1*age' or \n \
+ 'model=1+age+V1+V2' or 'model=1+age+V1+V2+V1*V2' etc. \n");		\
     if(mle != -1){
-      printf("Fix the model line and run imach with mle=-1 to get a correct template of the parameter file.\n");
+      printf("Fix the model line and run imach with mle=-1 to get a correct template of the parameter vectors and subdiagonal covariance matrix.\n");
       exit(1);
     }
   }
@@ -11145,16 +11300,6 @@ Please run with mle=-1 to get a correct
     
     fflush(ficlog);
     
-    /*-------- Rewriting parameter file ----------*/
-    strcpy(rfileres,"r");    /* "Rparameterfile */
-    strcat(rfileres,optionfilefiname);    /* Parameter file first name*/
-    strcat(rfileres,".");    /* */
-    strcat(rfileres,optionfilext);    /* Other files have txt extension */
-    if((ficres =fopen(rfileres,"w"))==NULL) {
-      printf("Problem writing new parameter file: %s\n", rfileres);goto end;
-      fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end;
-    }
-    fprintf(ficres,"#%s\n",version);
   }    /* End of mle != -3 */
   
   /*  Main data
@@ -11788,7 +11933,7 @@ Please run with mle=-1 to get a correct
     printf("\n");
     
     /*--------- results files --------------*/
-    fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, weightopt,model);
+    /* fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, weightopt,model); */
     
     
     fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
@@ -12586,6 +12731,8 @@ Please run with mle=-1 to get a correct
   fclose(ficlog);
   /*------ End -----------*/
   
+
+/* Executes gnuplot */
   
   printf("Before Current directory %s!\n",pathcd);
 #ifdef WIN32
@@ -12654,4 +12801,6 @@ end:
     printf("\nType  q for exiting: "); fflush(stdout);
     scanf("%s",z);
   }
+  printf("End\n");
+  exit(0);
 }