]> henry.ined.fr Git - .git/commitdiff
Summary: Better
authorN. Brouard <brouard@ined.fr>
Mon, 29 Aug 2016 07:53:18 +0000 (07:53 +0000)
committerN. Brouard <brouard@ined.fr>
Mon, 29 Aug 2016 07:53:18 +0000 (07:53 +0000)
src/imach.c

index ce6838fa3725bef5a8b5208ccdde95c848b0003e..00b0480bcd201f103f9911757ba0372bd5816c5c 100644 (file)
@@ -1,6 +1,11 @@
 /* $Id$
   $State$
   $Log$
+  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
 
@@ -1032,7 +1037,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 
 
@@ -2126,16 +2132,13 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
        if (k != i) {
          printf("%d%d ",i,k);
          fprintf(ficlog,"%d%d ",i,k);
-         fprintf(ficres,"%1d%1d ",i,k);
          for(j=1; j <=ncovmodel; j++){
            printf("%12.7f ",p[jk]);
            fprintf(ficlog,"%12.7f ",p[jk]);
-           fprintf(ficres,"%12.7f ",p[jk]);
            jk++; 
          }
          printf("\n");
          fprintf(ficlog,"\n");
-         fprintf(ficres,"\n");
        }
       }
     }
@@ -2256,7 +2259,7 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
       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]; 
@@ -4130,7 +4133,7 @@ void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **ag
     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));
@@ -4140,55 +4143,54 @@ void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **ag
   }
   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;
     }
@@ -4198,7 +4200,7 @@ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age
     /*         meanqt[m][z1]=0.; */
     /*   } */
     /* } */
-               
+    
     dateintsum=0;
     k2cpt=0;
     /* For that combination of covariate j1, we count and print the frequencies in one pass */
@@ -4277,35 +4279,41 @@ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age
     } /* 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> ");
@@ -4316,115 +4324,115 @@ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age
       }
     }
     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");
@@ -4442,7 +4450,7 @@ Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age
     fprintf(ficresphtmfr,"</table>\n");
   } /* end selected combination of covariate j1 */
   dateintmean=dateintsum/k2cpt; 
-       
+  
   fclose(ficresp);
   fclose(ficresphtm);
   fclose(ficresphtmfr);
@@ -6350,10 +6358,7 @@ void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar,
       
        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,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
       
        for (i=1; i<= nlstate ; i ++) {
          if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
@@ -6467,7 +6472,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(file
   /*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 ++) {
@@ -7827,14 +7832,36 @@ int readdata(char datafile[], int firstobs, int lastobs, int *imax)
   /*-------- 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);
@@ -8455,27 +8482,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy
 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;
@@ -8500,7 +8507,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy
       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++;
@@ -8553,167 +8560,167 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy
       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]);
@@ -10045,8 +10052,6 @@ Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numline
   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.
@@ -10864,10 +10869,12 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
       }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 */
@@ -10897,7 +10904,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
        break;
       else{ /* Processess output results for this combination of covariate values */
       }                                   
-    }
+    } /* end while */