--- imach/src/imach.c	2016/08/26 14:23:35	1.238
+++ imach/src/imach.c	2016/08/29 17:17:25	1.241
@@ -1,6 +1,17 @@
-/* $Id: imach.c,v 1.238 2016/08/26 14:23:35 brouard Exp $
+/* $Id: imach.c,v 1.241 2016/08/29 17:17:25 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.241  2016/08/29 17:17:25  brouard
+  Summary: gnuplot problem in Back projection to fix
+
+  Revision 1.240  2016/08/29 07:53:18  brouard
+  Summary: Better
+
+  Revision 1.239  2016/08/26 15:51:03  brouard
+  Summary: Improvement in Powell output in order to copy and paste
+
+  Author:
+
   Revision 1.238  2016/08/26 14:23:35  brouard
   Summary: Starting tests of 0.99
 
@@ -914,12 +925,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.238 2016/08/26 14:23:35 brouard Exp $ */
+/* $Id: imach.c,v 1.241 2016/08/29 17:17:25 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.238 $ $Date: 2016/08/26 14:23:35 $"; 
+char fullversion[]="$Revision: 1.241 $ $Date: 2016/08/29 17:17:25 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1032,7 +1043,8 @@ double dval;
 #define FTOL 1.0e-10
 
 #define NRANSI 
-#define ITMAX 200 
+#define ITMAX 200
+#define ITPOWMAX 20 /* This is now multiplied by the number of parameters */ 
 
 #define TOL 2.0e-4 
 
@@ -2061,10 +2073,10 @@ void powell(double p[], double **xi, int
  void linmin(double p[], double xi[], int n, double *fret, 
 	      double (*func)(double [])); 
 #else 
- void linmin(double p[], double xi[], int n, double *fret, 
-						 double (*func)(double []),int *flat); 
+ void linmin(double p[], double xi[], int n, double *fret,
+	     double (*func)(double []),int *flat); 
 #endif
-  int i,ibig,j; 
+ int i,ibig,j,jk,k; 
   double del,t,*pt,*ptt,*xit;
   double directest;
   double fp,fptt;
@@ -2096,31 +2108,64 @@ void powell(double p[], double **xi, int
     fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
 /*     fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */
     for (i=1;i<=n;i++) {
-      printf(" %d %.12f",i, p[i]);
-      fprintf(ficlog," %d %.12lf",i, p[i]);
       fprintf(ficrespow," %.12lf", p[i]);
     }
+    fprintf(ficrespow,"\n");fflush(ficrespow);
+    printf("\n#model=  1      +     age ");
+    fprintf(ficlog,"\n#model=  1      +     age ");
+    if(nagesqr==1){
+	printf("  + age*age  ");
+	fprintf(ficlog,"  + age*age  ");
+    }
+    for(j=1;j <=ncovmodel-2;j++){
+      if(Typevar[j]==0) {
+	printf("  +      V%d  ",Tvar[j]);
+	fprintf(ficlog,"  +      V%d  ",Tvar[j]);
+      }else if(Typevar[j]==1) {
+	printf("  +    V%d*age ",Tvar[j]);
+	fprintf(ficlog,"  +    V%d*age ",Tvar[j]);
+      }else if(Typevar[j]==2) {
+	printf("  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+	fprintf(ficlog,"  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
+      }
+    }
     printf("\n");
+/*     printf("12   47.0114589    0.0154322   33.2424412    0.3279905    2.3731903  */
+/* 13  -21.5392400    0.1118147    1.2680506    1.2973408   -1.0663662  */
     fprintf(ficlog,"\n");
-    fprintf(ficrespow,"\n");fflush(ficrespow);
-    if(*iter <=3){
+    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(j=1; j <=ncovmodel; j++){
+	    printf("%12.7f ",p[jk]);
+	    fprintf(ficlog,"%12.7f ",p[jk]);
+	    jk++; 
+	  }
+	  printf("\n");
+	  fprintf(ficlog,"\n");
+	}
+      }
+    }
+    if(*iter <=3 && *iter >1){
       tml = *localtime(&rcurr_time);
       strcpy(strcurr,asctime(&tml));
       rforecast_time=rcurr_time; 
       itmp = strlen(strcurr);
       if(strcurr[itmp-1]=='\n')  /* Windows outputs with a new line */
-				strcurr[itmp-1]='\0';
+	strcurr[itmp-1]='\0';
       printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
       fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
       for(niterf=10;niterf<=30;niterf+=10){
-				rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);
-				forecast_time = *localtime(&rforecast_time);
-				strcpy(strfor,asctime(&forecast_time));
-				itmp = strlen(strfor);
-				if(strfor[itmp-1]=='\n')
-					strfor[itmp-1]='\0';
-				printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
-				fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
+	rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);
+	forecast_time = *localtime(&rforecast_time);
+	strcpy(strfor,asctime(&forecast_time));
+	itmp = strlen(strfor);
+	if(strfor[itmp-1]=='\n')
+	  strfor[itmp-1]='\0';
+	printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
+	fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
       }
     }
     for (i=1;i<=n;i++) { /* For each direction i */
@@ -2220,7 +2265,7 @@ void powell(double p[], double **xi, int
       free_vector(pt,1,n); 
       return; 
     } /* enough precision */ 
-    if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); 
+    if (*iter == ITMAX*n) nrerror("powell exceeding maximum iterations."); 
     for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */
       ptt[j]=2.0*p[j]-pt[j]; 
       xit[j]=p[j]-pt[j]; 
@@ -4094,7 +4139,7 @@ void  freqsummary(char fileres[], int ia
     fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
     exit(0);
   }
-
+  
   strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));
   if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {
     printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
@@ -4104,55 +4149,54 @@ void  freqsummary(char fileres[], int ia
   }
   else{
     fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
-<hr size=\"2\" color=\"#EC5E5E\"> \n\
+<hr size=\"2\" color=\"#EC5E5E\"> \n					\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
 	    fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }
   fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm);
-    
+  
   strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
   if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
     printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
     fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
     fflush(ficlog);
     exit(70); 
-  }
-  else{
+  } else{
     fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
-<hr size=\"2\" color=\"#EC5E5E\"> \n\
+<hr size=\"2\" color=\"#EC5E5E\"> \n					\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
 	    fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }
-  fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
-
+  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);
   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
      med_educ                V1=1 V2=0, 
      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++)  
-				for(m=iagemin; m <= iagemax+3; m++)
-					freq[i][jk][m]=0;
-		
+	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;
+	prop[i][m]=0;
       posprop[i]=0;
       pospropt[i]=0;
     }
@@ -4162,7 +4206,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
     /* 	meanqt[m][z1]=0.; */
     /*   } */
     /* } */
-		
+    
     dateintsum=0;
     k2cpt=0;
     /* For that combination of covariate j1, we count and print the frequencies in one pass */
@@ -4241,35 +4285,41 @@ Title=%s <br>Datafile=%s Firstpass=%d La
     } /* 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  (ncoveff>0) { */
-    if  (cptcoveff>0) {
+    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++){
-	fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-	fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
-	fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+	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#********** Variable "); 
-      for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
       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(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> ");
@@ -4280,115 +4330,115 @@ Title=%s <br>Datafile=%s Firstpass=%d La
       }
     }
     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> ");
+	fprintf(ficlog,"1");
+	fprintf(ficresphtmfr,"<tr><th>0</th> ");
       }else if(iage==iagemax+2){
-				fprintf(ficlog,"0");
-				fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
+	fprintf(ficlog,"0");
+	fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
       }else if(iage==iagemax+3){
-				fprintf(ficlog,"Total");
-				fprintf(ficresphtmfr,"<tr><th>Total</th> ");
+	fprintf(ficlog,"Total");
+	fprintf(ficresphtmfr,"<tr><th>Total</th> ");
       }else{
-				if(first==1){
-					first=0;
-					printf("See log file for details...\n");
-				}
-				fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage);
-				fprintf(ficlog,"Age %d", iage);
+	if(first==1){
+	  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(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);
-				}
+	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);
+	}
       }
-			
+      
       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];
+	/* 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] */
+	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] */
       }
       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){
-					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);
-					}
-				}
-				pospropt[jk] +=posprop[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){
+	  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);
+	  }
+	}
+	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]);
-						}
-						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]);
-				}
+	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]);
+	    }
+	    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];
+	posproptt += pospropt[jk];
       }
       fprintf(ficresphtmfr,"</tr>\n ");
       if(iage <= iagemax){
-				fprintf(ficresp,"\n");
-				fprintf(ficresphtm,"</tr>\n");
+	fprintf(ficresp,"\n");
+	fprintf(ficresphtm,"</tr>\n");
       }
       if(first==1)
-				printf("Others in log...\n");
+	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);	
+	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,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);	
       }
     }
     fprintf(ficresphtm,"</tr>\n");
@@ -4406,7 +4456,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
     fprintf(ficresphtmfr,"</table>\n");
   } /* end selected combination of covariate j1 */
   dateintmean=dateintsum/k2cpt; 
-	
+  
   fclose(ficresp);
   fclose(ficresphtm);
   fclose(ficresphtmfr);
@@ -5568,7 +5618,9 @@ void  concatwav(int wav[], int **dh, int
   
   pstamp(ficresvpl);
   fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n");
-  fprintf(ficresvpl,"# Age");
+  fprintf(ficresvpl,"# Age ");
+  if(nresult >=1)
+    fprintf(ficresvpl," Result# ");
   for(i=1; i<=nlstate;i++)
       fprintf(ficresvpl," %1d-%1d",i,i);
   fprintf(ficresvpl,"\n");
@@ -5654,6 +5706,8 @@ void  concatwav(int wav[], int **dh, int
       varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */
 
     fprintf(ficresvpl,"%.0f ",age );
+    if(nresult >=1)
+      fprintf(ficresvpl,"%d ",nres );
     for(i=1; i<=nlstate;i++)
       fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age]));
     fprintf(ficresvpl,"\n");
@@ -6058,7 +6112,7 @@ void printinghtml(char fileresu[], char
    jj1=0;
 
    for(nres=1; nres <= nresult; nres++) /* For each resultline */
-   for(k1=1; k1<=m;k1++){
+   for(k1=1; k1<=m;k1++){ /* For each combination of covariate */
      if(TKresult[nres]!= k1)
        continue;
 
@@ -6086,51 +6140,51 @@ 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.svg\">%s_%d-1.svg</a><br> \
-<img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
+     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> \
+<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.svg\">%s_%d-2.svg</a><br> \
-<img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);     
+     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> \
+<img src=\"%s_%d-2-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);     
      /* Quasi-incidences */
      fprintf(fichtm,"<br>\n- I<sub>ij</sub> or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
  before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too, \
  incidence (rates) are the limit when h tends to zero of the ratio of the probability  <sub>h</sub>P<sub>ij</sub> \
-divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \
-<img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); 
+divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3-%d.svg\">%s_%d-3-%d.svg</a><br> \
+<img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); 
      /* Survival functions (period) in state j */
      for(cpt=1; cpt<=nlstate;cpt++){
-       fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
-<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1);
+       fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
+<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
      }
      /* State specific survival functions (period) */
      for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\
  Or probability to survive in various states (1 to %d) being in state %d at different ages.	\
- <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> <img src=\"%s_%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1);
+ <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++){
-       fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
-<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);
+       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> \
+<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.svg\">%s_%d-%d.svg</a><br> \
-<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1);
+	 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> \
+<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.svg\">%s%d_%d.svg</a><br> \
-<img src=\"%s_%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1);
+	 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);
        }
      }
 	 
      for(cpt=1; cpt<=nlstate;cpt++) {
-       fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \
-<img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1);
+       fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a> <br> \
+<img src=\"%s_%d-%d-%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres);
      }
      /* } /\* end i1 *\/ */
    }/* End k1 */
@@ -6189,7 +6243,7 @@ See page 'Matrix of variance-covariance
 
    jj1=0;
 
-   for(nres=1; nres <= nresult; nres++) /* For each resultline */
+   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
    for(k1=1; k1<=m;k1++){
      if(TKresult[nres]!= k1)
        continue;
@@ -6213,17 +6267,18 @@ See page 'Matrix of variance-covariance
      }
      for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \
-prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d.svg\"> %s_%d-%d.svg</a>\n <br>\
-<img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1);  
+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);  
      }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and \
 health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
 true period expectancies (those weighted with period prevalences are also\
  drawn in addition to the population based expectancies computed using\
- observed and cahotic prevalences:  <a href=\"%s_%d.svg\">%s_%d.svg</a>\n<br>\
-<img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1);
+ observed and cahotic prevalences:  <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>\n<br>\
+<img src=\"%s_%d-%d.svg\">",subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres);
      /* } /\* end i1 *\/ */
    }/* End k1 */
+  }/* End nres */
    fprintf(fichtm,"</ul>");
    fflush(fichtm);
 }
@@ -6312,23 +6367,20 @@ void printinggnuplot(char fileresu[], ch
 	  continue;
 	}
       
-	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
-	fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
-	fprintf(ficgp,"set xlabel \"Age\" \n\
-set ylabel \"Probability\" \n		  \
-set ter svg size 640, 480\n						\
-plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+	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);
       
 	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+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+	fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $4+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-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-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); 
+	fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2==%d ? $4-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); 
 	for (i=1; i<= nlstate ; i ++) {
 	  if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
 	  else fprintf(ficgp," %%*lf (%%*lf)");
@@ -6336,7 +6388,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 	fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
 	if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */
 	  /* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */
-	  fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */
+	  fprintf(ficgp,",\"%s\" u ($2==%d ?$1:1/0):(",subdirf2(fileresu,"PLB_"),nres); /* Age is in 1, nres in 2 to be fixed */
 	  if(cptcoveff ==0){
 	    fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ",	 2+(cpt-1),  cpt );
 	  }else{
@@ -6393,7 +6445,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 	continue;
       }
 			
-      fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
+      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1,nres);
       for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
 	if(vpopbased==0)
 	  fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
@@ -6431,7 +6483,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
   /*3eme*/
   for (k1=1; k1<= m ; k1 ++){
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-      if(TKresult[nres]!= k)
+      if(TKresult[nres]!= k1)
 	continue;
 
       for (cpt=1; cpt<= nlstate ; cpt ++) {
@@ -6455,7 +6507,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 			
 	/*       k=2+nlstate*(2*cpt-2); */
 	k=2+(nlstate+1)*(cpt-1);
-	fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
+	fprintf(ficgp,"\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);
 	/*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
@@ -6501,7 +6553,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 	  continue;
 	}
       
-	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
+	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
 	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
 set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
 	k=3;
@@ -6547,7 +6599,7 @@ set ter svg size 640, 480\nunset log y\n
 	  continue;
 	}
       
-	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
+	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
 	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
 set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
 	k=3;
@@ -6602,7 +6654,7 @@ set ter svg size 640, 480\nunset log y\n
 	continue;
       }
       
-      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);
+      fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
       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 */
@@ -6648,7 +6700,7 @@ set ter svg size 640, 480\nunset log y\n
 	  continue;
 	}
 	
-	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
+	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
 	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 */
@@ -6700,7 +6752,7 @@ set ter svg size 640, 480\nunset log y\n
 	}
 	
 	fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
-	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
+	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
 	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
 set ter svg size 640, 480\nunset log y\nplot [%.f:%.f]  ", ageminpar, agemaxpar);
 	for (i=1; i<= nlstate+1 ; i ++){  /* nlstate +1 p11 p21 p.1 */
@@ -6816,7 +6868,7 @@ set ter svg size 640, 480\nunset log y\n
 	fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
       }	
       fprintf(ficgp,"\n#\n");
-      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng);
+      fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng,nres);
       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 */
@@ -7791,14 +7843,36 @@ int readdata(char datafile[], int firsto
   /*-------- data file ----------*/
   FILE *fic;
   char dummy[]="                         ";
-  int i=0, j=0, n=0, iv=0;
+  int i=0, j=0, n=0, iv=0, v;
   int lstra;
   int linei, month, year,iout;
   char line[MAXLINE], linetmp[MAXLINE];
   char stra[MAXLINE], strb[MAXLINE];
   char *stratrunc;
 
+  DummyV=ivector(1,NCOVMAX); /* 1 to 3 */
+  FixedV=ivector(1,NCOVMAX); /* 1 to 3 */
 
+  for(v=1; v <=ncovcol;v++){
+    DummyV[v]=0;
+    FixedV[v]=0;
+  }
+  for(v=ncovcol+1; v <=ncovcol+nqv;v++){
+    DummyV[v]=1;
+    FixedV[v]=0;
+  }
+  for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){
+    DummyV[v]=0;
+    FixedV[v]=1;
+  }
+  for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){
+    DummyV[v]=1;
+    FixedV[v]=1;
+  }
+  for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){
+    printf("Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
+    fprintf(ficlog,"Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
+  }
 
   if((fic=fopen(datafile,"r"))==NULL)    {
     printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);
@@ -8419,27 +8493,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(v=1; v <=ncovcol;v++){
-    DummyV[v]=0;
-    FixedV[v]=0;
-  }
-  for(v=ncovcol+1; v <=ncovcol+nqv;v++){
-    DummyV[v]=1;
-    FixedV[v]=0;
-  }
-  for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){
-    DummyV[v]=0;
-    FixedV[v]=1;
-  }
-  for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){
-    DummyV[v]=1;
-    FixedV[v]=1;
-  }
-  for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){
-    printf("Decodemodel: V%d, Dummy(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
-    fprintf(ficlog,"Decodemodel: V%d, Dummy(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
-  }
+  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;
@@ -8464,7 +8518,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       TvarFind[ncovf]=k;
       TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
       TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
-    }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/ /* Only simple fixed quantitative variable */
+    }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){/* Remind that product Vn*Vm are added in k Only simple fixed quantitative variable */
       Fixed[k]= 0;
       Dummy[k]= 1;
       nqfveff++;
@@ -8517,167 +8571,167 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
       TvarA[ncova]=Tvar[k];
       TvarAind[ncova]=k;
       if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */
-	Fixed[k]= 2;
-	Dummy[k]= 2;
-	modell[k].maintype= ATYPE;
-	modell[k].subtype= APFD;
-	/* ncoveff++; */
+      	Fixed[k]= 2;
+      	Dummy[k]= 2;
+      	modell[k].maintype= ATYPE;
+      	modell[k].subtype= APFD;
+      	/* ncoveff++; */
       }else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product Vn*Vm are added in k*/
-	Fixed[k]= 2;
-	Dummy[k]= 3;
-	modell[k].maintype= ATYPE;
-	modell[k].subtype= APFQ;		/*	Product age * fixed quantitative */
-	/* nqfveff++;  /\* Only simple fixed quantitative variable *\/ */
+      	Fixed[k]= 2;
+      	Dummy[k]= 3;
+      	modell[k].maintype= ATYPE;
+      	modell[k].subtype= APFQ;		/*	Product age * fixed quantitative */
+      	/* nqfveff++;  /\* Only simple fixed quantitative variable *\/ */
       }else if( Tvar[k] <=ncovcol+nqv+ntv ){
-	Fixed[k]= 3;
-	Dummy[k]= 2;
-	modell[k].maintype= ATYPE;
-	modell[k].subtype= APVD;		/*	Product age * varying dummy */
-	/* ntveff++; /\* Only simple time varying dummy variable *\/ */
+      	Fixed[k]= 3;
+      	Dummy[k]= 2;
+      	modell[k].maintype= ATYPE;
+      	modell[k].subtype= APVD;		/*	Product age * varying dummy */
+      	/* ntveff++; /\* Only simple time varying dummy variable *\/ */
       }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){
-	Fixed[k]= 3;
-	Dummy[k]= 3;
-	modell[k].maintype= ATYPE;
-	modell[k].subtype= APVQ;		/*	Product age * varying quantitative */
-	/* nqtveff++;/\* Only simple time varying quantitative variable *\/ */
+      	Fixed[k]= 3;
+      	Dummy[k]= 3;
+      	modell[k].maintype= ATYPE;
+      	modell[k].subtype= APVQ;		/*	Product age * varying quantitative */
+      	/* nqtveff++;/\* Only simple time varying quantitative variable *\/ */
       }
     }else if (Typevar[k] == 2) {  /* product without age */
       k1=Tposprod[k];
       if(Tvard[k1][1] <=ncovcol){
-	if(Tvard[k1][2] <=ncovcol){
-	  Fixed[k]= 1;
-	  Dummy[k]= 0;
-	  modell[k].maintype= FTYPE;
-	  modell[k].subtype= FPDD;		/*	Product fixed dummy * fixed dummy */
-	  ncovf++; /* Fixed variables without age */
-	  TvarF[ncovf]=Tvar[k];
-	  TvarFind[ncovf]=k;
-	}else if(Tvard[k1][2] <=ncovcol+nqv){
-	  Fixed[k]= 0;  /* or 2 ?*/
-	  Dummy[k]= 1;
-	  modell[k].maintype= FTYPE;
-	  modell[k].subtype= FPDQ;		/*	Product fixed dummy * fixed quantitative */
-	  ncovf++; /* Varying variables without age */
-	  TvarF[ncovf]=Tvar[k];
-	  TvarFind[ncovf]=k;
-	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
-	  Fixed[k]= 1;
-	  Dummy[k]= 0;
-	  modell[k].maintype= VTYPE;
-	  modell[k].subtype= VPDD;		/*	Product fixed dummy * varying dummy */
-	  ncovv++; /* Varying variables without age */
-	  TvarV[ncovv]=Tvar[k];
-	  TvarVind[ncovv]=k;
-	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
-	  Fixed[k]= 1;
-	  Dummy[k]= 1;
-	  modell[k].maintype= VTYPE;
-	  modell[k].subtype= VPDQ;		/*	Product fixed dummy * varying quantitative */
-	  ncovv++; /* Varying variables without age */
-	  TvarV[ncovv]=Tvar[k];
-	  TvarVind[ncovv]=k;
-	} 
+    	if(Tvard[k1][2] <=ncovcol){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 0;
+    	  modell[k].maintype= FTYPE;
+    	  modell[k].subtype= FPDD;		/*	Product fixed dummy * fixed dummy */
+    	  ncovf++; /* Fixed variables without age */
+    	  TvarF[ncovf]=Tvar[k];
+    	  TvarFind[ncovf]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv){
+    	  Fixed[k]= 0;  /* or 2 ?*/
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= FTYPE;
+    	  modell[k].subtype= FPDQ;		/*	Product fixed dummy * fixed quantitative */
+    	  ncovf++; /* Varying variables without age */
+    	  TvarF[ncovf]=Tvar[k];
+    	  TvarFind[ncovf]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 0;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDD;		/*	Product fixed dummy * varying dummy */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDQ;		/*	Product fixed dummy * varying quantitative */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}
       }else if(Tvard[k1][1] <=ncovcol+nqv){
-	if(Tvard[k1][2] <=ncovcol){
-	  Fixed[k]= 0;  /* or 2 ?*/
-	  Dummy[k]= 1;
-	  modell[k].maintype= FTYPE;
-	  modell[k].subtype= FPDQ;		/*	Product fixed quantitative * fixed dummy */
-	  ncovf++; /* Fixed variables without age */
-	  TvarF[ncovf]=Tvar[k];
-	  TvarFind[ncovf]=k;
-	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
-	  Fixed[k]= 1;
-	  Dummy[k]= 1;
-	  modell[k].maintype= VTYPE;
-	  modell[k].subtype= VPDQ;		/*	Product fixed quantitative * varying dummy */
-	  ncovv++; /* Varying variables without age */
-	  TvarV[ncovv]=Tvar[k];
-	  TvarVind[ncovv]=k;
-	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
-	  Fixed[k]= 1;
-	  Dummy[k]= 1;
-	  modell[k].maintype= VTYPE;
-	  modell[k].subtype= VPQQ;		/*	Product fixed quantitative * varying quantitative */
-	  ncovv++; /* Varying variables without age */
-	  TvarV[ncovv]=Tvar[k];
-	  TvarVind[ncovv]=k;
-	  ncovv++; /* Varying variables without age */
-	  TvarV[ncovv]=Tvar[k];
-	  TvarVind[ncovv]=k;
-	} 
+    	if(Tvard[k1][2] <=ncovcol){
+    	  Fixed[k]= 0;  /* or 2 ?*/
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= FTYPE;
+    	  modell[k].subtype= FPDQ;		/*	Product fixed quantitative * fixed dummy */
+    	  ncovf++; /* Fixed variables without age */
+    	  TvarF[ncovf]=Tvar[k];
+    	  TvarFind[ncovf]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDQ;		/*	Product fixed quantitative * varying dummy */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPQQ;		/*	Product fixed quantitative * varying quantitative */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}
       }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){
-	if(Tvard[k1][2] <=ncovcol){
-	  Fixed[k]= 1;
-	  Dummy[k]= 1;
-	  modell[k].maintype= VTYPE;
-	  modell[k].subtype= VPDD;		/*	Product time varying dummy * fixed dummy */
-	  ncovv++; /* Varying variables without age */
-	  TvarV[ncovv]=Tvar[k];
-	  TvarVind[ncovv]=k;
-	}else if(Tvard[k1][2] <=ncovcol+nqv){
-	  Fixed[k]= 1;
-	  Dummy[k]= 1;
-	  modell[k].maintype= VTYPE;
-	  modell[k].subtype= VPDQ;		/*	Product time varying dummy * fixed quantitative */
-	  ncovv++; /* Varying variables without age */
-	  TvarV[ncovv]=Tvar[k];
-	  TvarVind[ncovv]=k;
-	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
-	  Fixed[k]= 1;
-	  Dummy[k]= 0;
-	  modell[k].maintype= VTYPE;
-	  modell[k].subtype= VPDD;		/*	Product time varying dummy * time varying dummy */
-	  ncovv++; /* Varying variables without age */
-	  TvarV[ncovv]=Tvar[k];
-	  TvarVind[ncovv]=k;
-	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
-	  Fixed[k]= 1;
-	  Dummy[k]= 1;
-	  modell[k].maintype= VTYPE;
-	  modell[k].subtype= VPDQ;		/*	Product time varying dummy * time varying quantitative */
-	  ncovv++; /* Varying variables without age */
-	  TvarV[ncovv]=Tvar[k];
-	  TvarVind[ncovv]=k;
-	} 
+    	if(Tvard[k1][2] <=ncovcol){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDD;		/*	Product time varying dummy * fixed dummy */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDQ;		/*	Product time varying dummy * fixed quantitative */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 0;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDD;		/*	Product time varying dummy * time varying dummy */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDQ;		/*	Product time varying dummy * time varying quantitative */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}
       }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){
-	if(Tvard[k1][2] <=ncovcol){
-	  Fixed[k]= 1;
-	  Dummy[k]= 1;
-	  modell[k].maintype= VTYPE;
-	  modell[k].subtype= VPDQ;		/*	Product time varying quantitative * fixed dummy */
-	  ncovv++; /* Varying variables without age */
-	  TvarV[ncovv]=Tvar[k];
-	  TvarVind[ncovv]=k;
-	}else if(Tvard[k1][2] <=ncovcol+nqv){
-	  Fixed[k]= 1;
-	  Dummy[k]= 1;
-	  modell[k].maintype= VTYPE;
-	  modell[k].subtype= VPQQ;		/*	Product time varying quantitative * fixed quantitative */
-	  ncovv++; /* Varying variables without age */
-	  TvarV[ncovv]=Tvar[k];
-	  TvarVind[ncovv]=k;
-	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
-	  Fixed[k]= 1;
-	  Dummy[k]= 1;
-	  modell[k].maintype= VTYPE;
-	  modell[k].subtype= VPDQ;		/*	Product time varying quantitative * time varying dummy */
-	  ncovv++; /* Varying variables without age */
-	  TvarV[ncovv]=Tvar[k];
-	  TvarVind[ncovv]=k;
-	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
-	  Fixed[k]= 1;
-	  Dummy[k]= 1;
-	  modell[k].maintype= VTYPE;
-	  modell[k].subtype= VPQQ;		/*	Product time varying quantitative * time varying quantitative */
-	  ncovv++; /* Varying variables without age */
-	  TvarV[ncovv]=Tvar[k];
-	  TvarVind[ncovv]=k;
-	} 
+    	if(Tvard[k1][2] <=ncovcol){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDQ;		/*	Product time varying quantitative * fixed dummy */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPQQ;		/*	Product time varying quantitative * fixed quantitative */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPDQ;		/*	Product time varying quantitative * time varying dummy */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+    	  Fixed[k]= 1;
+    	  Dummy[k]= 1;
+    	  modell[k].maintype= VTYPE;
+    	  modell[k].subtype= VPQQ;		/*	Product time varying quantitative * time varying quantitative */
+    	  ncovv++; /* Varying variables without age */
+    	  TvarV[ncovv]=Tvar[k];
+    	  TvarVind[ncovv]=k;
+    	}
       }else{
-	printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
-	fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
-      } /* end k1 */
+    	printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
+    	fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
+      } /*end k1*/
     }else{
       printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
       fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
@@ -10009,8 +10063,6 @@ Please run with mle=-1 to get a correct
   Typevar=ivector(-1,NCOVMAX); /* -1 to 2 */
   Fixed=ivector(-1,NCOVMAX); /* -1 to 3 */
   Dummy=ivector(-1,NCOVMAX); /* -1 to 3 */
-  DummyV=ivector(1,NCOVMAX); /* 1 to 3 */
-  FixedV=ivector(1,NCOVMAX); /* 1 to 3 */
   /*  V2+V1+V4+age*V3 is a model with 4 covariates (3 plus signs). 
       For each model-covariate stores the data-covariate id. Tvar[1]=2, Tvar[2]=1, Tvar[3]=4, 
       Tvar[4=age*V3] is 3 and 'age' is recorded in Tage.
@@ -10828,10 +10880,12 @@ Please run with mle=-1 to get a correct
       }else
 	break;
     }
+    if (!feof(ficpar))
     while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){
-      if (num_filled == 0)
+      if (num_filled == 0){
 	resultline[0]='\0';
-      else if (num_filled != 1){
+      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 */
@@ -10861,7 +10915,7 @@ Please run with mle=-1 to get a correct
 	break;
       else{ /* Processess output results for this combination of covariate values */
       }				   
-    }
+    } /* end while */