Summary: Big bug thanks to Flavia

Even model=1+age+V2. did not work anymore

diff --git a/src/imach.c b/src/imach.c
index c7fe66de87a87c1234617da288dcecb9de158ee2..56d2e1192d85e5e92aff08d27c0735874cccb589 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -1,6 +1,9 @@
 /* $Id$
   $State$
   $Log$
+  Revision 1.199  2015/09/07 14:09:23  brouard
+  Summary: 0.98q6 changing default small png format for graph to vectorized svg.
+
   Revision 1.198  2015/09/03 07:14:39  brouard
   Summary: 0.98q5 Flavia
 
@@ -878,7 +881,7 @@ double  idx;
 int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
 int *Tage;
 int *Ndum; /** Freq of modality (tricode */
-int **codtab; /**< codtab=imatrix(1,100,1,10); */
+/* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */
 int **Tvard, *Tprod, cptcovprod, *Tvaraff;
 double *lsurv, *lpop, *tpop;
 
@@ -1914,13 +1917,16 @@ double **prevalim(double **prlim, int nlstate, double x[], double age, double **
     if(nagesqr==1)
       cov[3]= agefin*agefin;;
     for (k=1; k<=cptcovn;k++) {
-      cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])];
+      /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+      cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
       /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */
     }
     /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-    for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2];
+    /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */
+    for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2];
     for (k=1; k<=cptcovprod;k++) /* Useless */
-      cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])];
+      /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+      cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
     
     /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
     /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
@@ -2093,12 +2099,15 @@ double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, in
       if(nagesqr==1)
        cov[3]= agexact*agexact;
       for (k=1; k<=cptcovn;k++) 
-       cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])];
+       cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+       /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
       for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
        /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-       cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2];
+       cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+       /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
       for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
-       cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])];
+       cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
+       /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
 
 
       /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
@@ -3839,7 +3848,8 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
   }  
   fprintf(ficresprobmorprev,"\n");
   fprintf(ficgp,"\n# Routine varevsij");
-  /* fprintf(fichtm, "#Local time at start: %s", strstart);*/
+  fprintf(ficgp,"\nunset title \n");
+/* fprintf(fichtm, "#Local time at start: %s", strstart);*/
   fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
   fprintf(fichtm,"\n<br>%s  <br>\n",digitp);
 /*   } */
@@ -4055,6 +4065,7 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
   fprintf(ficgp,"\nunset parametric;unset label; set ter svg size 640, 480");
   /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
   fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
+  fprintf(ficgp,"\nset out \"%s%s.svg\";",subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
 /*   fprintf(ficgp,"\n plot \"%s\"  u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */
 /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
 /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
@@ -4066,7 +4077,7 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
   /*  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.svg\"> <br>\n", stepm,YEARM,digitp,digit);
 */
 /*   fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.svg\";replot;",digitp,optionfilefiname,digit); */
-  fprintf(ficgp,"\nset out \"%s%s.svg\";replot;\n",subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
+  fprintf(ficgp,"\nset out;\nset out \"%s%s.svg\";replot;set out;\n",subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
 
   free_vector(xp,1,npar);
   free_matrix(doldm,1,nlstate,1,nlstate);
@@ -4239,7 +4250,7 @@ void varprob(char optionfilefiname[], double **matcov, double x[], double delti[
   fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
   fprintf(fichtm,"\n");
 
-  fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4></br>this page is important in order to visualize confidence intervals and especially correlation between disability and recovery</li>\n",optionfilehtmcov);
+  fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.</li>\n",optionfilehtmcov);
   fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov);
   fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \
 and drawn. It helps understanding how is the covariance between two incidences.\
@@ -4291,7 +4302,8 @@ To be simple, these graphs help to understand the significativity of each parame
        if(nagesqr==1)
          cov[3]= age*age;
        for (k=1; k<=cptcovn;k++) {
-         cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];/* j1 1 2 3 4
+         cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)];
+         /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4
                                                         * 1  1 1 1 1
                                                         * 2  2 1 1 1
                                                         * 3  1 2 1 1
@@ -4299,9 +4311,9 @@ To be simple, these graphs help to understand the significativity of each parame
          /* nbcode[1][1]=0 nbcode[1][2]=1;*/
        }
        /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-       for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2];
+       for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
        for (k=1; k<=cptcovprod;k++)
-         cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])];
+         cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
        
     
        for(theta=1; theta <=npar; theta++){
@@ -4449,6 +4461,7 @@ To be simple, these graphs help to understand the significativity of each parame
                  /* mu2+ v21*lc1*cost + v22*lc2*sin(t) */
                  if(first==1){
                    first=0;
+                   fprintf(ficgp,"\n# Ellipsoids of confidence\n#\n");
                    fprintf(ficgp,"\nset parametric;unset label");
                    fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2);
                    fprintf(ficgp,"\nset ter svg size 640, 480");
@@ -4476,7 +4489,7 @@ To be simple, these graphs help to understand the significativity of each parame
                  }/* if first */
                } /* age mod 5 */
              } /* end loop age */
-             fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.svg\";replot;",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+             fprintf(ficgp,"\nset out;\nset out \"%s%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
              first=1;
            } /*l12 */
          } /* k12 */
@@ -4787,7 +4800,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar);
   fprintf(ficgp,"##############\n#\n");
 
   /*goto avoid;*/
-  fprintf(ficgp,"\n##############\n#Graphics of of probabilities or incidences\n#############\n");
+  fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n");
   fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n");
   fprintf(ficgp,"# logi(p12/p11)=p1 +p2*age +p3*age*age+ p4*V1+ p5*V1*age\n");
   fprintf(ficgp,"# logi(p13/p11)=a13+b13*age+c13age*age+d13*V1+e13*V1*age\n");
@@ -4810,7 +4823,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar);
        if (ng==2)
         fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
        else
-        fprintf(ficgp,"\nset title \"Probability\"\n");
+        fprintf(ficgp,"\nunset title \n");
        fprintf(ficgp,"\nset ter svg size 640, 480\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
        i=1;
        for(k2=1; k2<=nlstate; k2++) {
@@ -4832,7 +4845,8 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar);
               /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
               if(ij <=cptcovage) { /* Bug valgrind */
                 if((j-2)==Tage[ij]) { /* Bug valgrind */
-                  fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]);
+                  fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+                  /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
                   ij++;
                 }
               }
@@ -4851,7 +4865,8 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar);
               for(j=3; j <=ncovmodel-nagesqr; j++){
                 if(ij <=cptcovage) { /* Bug valgrind */
                   if((j-2)==Tage[ij]) { /* Bug valgrind */
-                    fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]);
+                    fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+                    /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
                     ij++;
                   }
                 }
@@ -6176,9 +6191,9 @@ int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxp
        k=k+1;
        /* to clean */
        //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
-       fprintf(ficrespl,"\n#******");
-       printf("\n#******");
-       fprintf(ficlog,"\n#******");
+       fprintf(ficrespl,"#******");
+       printf("#******");
+       fprintf(ficlog,"#******");
        for(j=1;j<=cptcoveff;j++) {
          fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
          printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
@@ -6190,7 +6205,7 @@ int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxp
 
        fprintf(ficrespl,"#Age ");
        for(j=1;j<=cptcoveff;j++) {
-         fprintf(ficrespl,"V%d %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+         fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
        }
        for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
        fprintf(ficrespl,"\n");
@@ -6918,7 +6933,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numline
      V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/
   /* 1 to ncodemax[j] is the maximum value of this jth covariate */
 
-  codtab=imatrix(1,100,1,10); /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
+  /*  codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
   /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/
   /* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/
   h=0;
@@ -6953,10 +6968,10 @@ Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numline
           */
   for(h=1; h <=100 ;h++){ 
     /* printf("h=%2d ", h); */
-     for(k=1; k <=10; k++){
+     /* for(k=1; k <=10; k++){ */
        /* printf("k=%d %d ",k,codtabm(h,k)); */
-       codtab[h][k]=codtabm(h,k);
-     }
+     /*   codtab[h][k]=codtabm(h,k); */
+     /* } */
      /* printf("\n"); */
   }
   /* for(k=1;k<=cptcoveff; k++){ /\* scans any effective covariate *\/ */
@@ -7650,7 +7665,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
     for (k=1; k <= (int) pow(2,cptcoveff); k++){
        fprintf(ficreseij,"\n#****** ");
        for(j=1;j<=cptcoveff;j++) {
-         fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+         fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
        }
        fprintf(ficreseij,"******\n");
 
@@ -7710,21 +7725,21 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
     for (k=1; k <= (int) pow(2,cptcoveff); k++){
        fprintf(ficrest,"\n#****** ");
        for(j=1;j<=cptcoveff;j++) 
-         fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+         fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
        fprintf(ficrest,"******\n");
 
        fprintf(ficresstdeij,"\n#****** ");
        fprintf(ficrescveij,"\n#****** ");
        for(j=1;j<=cptcoveff;j++) {
-         fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-         fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+         fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+         fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
        }
        fprintf(ficresstdeij,"******\n");
        fprintf(ficrescveij,"******\n");
 
        fprintf(ficresvij,"\n#****** ");
        for(j=1;j<=cptcoveff;j++) 
-         fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+         fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
        fprintf(ficresvij,"******\n");
 
        eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
@@ -7820,7 +7835,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
     for (k=1; k <= (int) pow(2,cptcoveff); k++){
        fprintf(ficresvpl,"\n#****** ");
        for(j=1;j<=cptcoveff;j++) 
-         fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+         fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
        fprintf(ficresvpl,"******\n");
       
        varpl=matrix(1,nlstate,(int) bage, (int) fage);
@@ -7857,7 +7872,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
     free_ivector(Tage,1,NCOVMAX);
 
     free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
-    free_imatrix(codtab,1,100,1,10);
+    /* free_imatrix(codtab,1,100,1,10); */
   fflush(fichtm);
   fflush(ficgp);