Some errors corrected in diagonalisation of covariance-matrice

diff --git a/src/imach.c b/src/imach.c
index 1ec326d8183ddf97d82dd4a0a419a135c7cef755..57287d605d42941985beadc7f08707d3eaf1cf6f 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -2154,7 +2154,7 @@ void varprob(char optionfilefiname[], double **matcov, double x[], double delti[
   int k2, l2, j1,  z1;\r
   int k=0,l, cptcode;\r
   int first=1, first1;\r
-  double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2;\r
+  double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp;\r
   double **dnewm,**doldm;\r
   double *xp;\r
   double *gp, *gm;\r
@@ -2378,15 +2378,15 @@ void varprob(char optionfilefiname[], double **matcov, double x[], double delti[
 \r
       /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/\r
       first1=1;\r
-      for (k1=1; k1<=(nlstate);k1++){\r
-       for (l1=1; l1<=(nlstate+ndeath);l1++){ \r
-         if(l1==k1) continue;\r
-         i=(k1-1)*(nlstate+ndeath)+l1;\r
-         for (k2=1; k2<=(nlstate);k2++){\r
-           for (l2=1; l2<=(nlstate+ndeath);l2++){ \r
-             if(l2==k2) continue;\r
-             j=(k2-1)*(nlstate+ndeath)+l2;\r
-             if(j<=i) continue;\r
+      for (k2=1; k2<=(nlstate);k2++){\r
+       for (l2=1; l2<=(nlstate+ndeath);l2++){ \r
+         if(l2==k2) continue;\r
+         j=(k2-1)*(nlstate+ndeath)+l2;\r
+         for (k1=1; k1<=(nlstate);k1++){\r
+           for (l1=1; l1<=(nlstate+ndeath);l1++){ \r
+             if(l1==k1) continue;\r
+             i=(k1-1)*(nlstate+ndeath)+l1;\r
+             if(i<=j) continue;\r
              for (age=bage; age<=fage; age ++){ \r
                if ((int)age %5==0){\r
                  v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM;\r
@@ -2394,55 +2394,49 @@ void varprob(char optionfilefiname[], double **matcov, double x[], double delti[
                  cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM;\r
                  mu1=mu[i][(int) age]/stepm*YEARM ;\r
                  mu2=mu[j][(int) age]/stepm*YEARM;\r
+                 c12=cv12/sqrt(v1*v2);\r
                  /* Computing eigen value of matrix of covariance */\r
-                 lc1=(v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));\r
-                 lc2=(v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12));\r
-                 if(first1==1){\r
-                   first1=0;\r
-                   printf("Var %.4e %.4e cov %.4e Eigen %.3e %.3e\nOthers in log...\n",v1,v2,cv12,lc1,lc2);\r
-                 }\r
-                 fprintf(ficlog,"Var %.4e %.4e cov %.4e Eigen %.3e %.3e\n",v1,v2,cv12,lc1,lc2);\r
+                 lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;\r
+                 lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;\r
                  /* Eigen vectors */\r
                  v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));\r
-                 v21=sqrt(1.-v11*v11);\r
+                 /*v21=sqrt(1.-v11*v11); *//* error */\r
+                 v21=(lc1-v1)/cv12*v11;\r
                  v12=-v21;\r
                  v22=v11;\r
+                 tnalp=v21/v11;\r
+                 if(first1==1){\r
+                   first1=0;\r
+                   printf("%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tang %.3f\nOthers in log...\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp);\r
+                 }\r
+                 fprintf(ficlog,"%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tan %.3f\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp);\r
                  /*printf(fignu*/\r
                  /* mu1+ v11*lc1*cost + v12*lc2*sin(t) */\r
-                 /* mu2+ v21*lc1*cost + v21*lc2*sin(t) */\r
+                 /* mu2+ v21*lc1*cost + v22*lc2*sin(t) */\r
                  if(first==1){\r
                    first=0;\r
-                   fprintf(ficgp,"\nset parametric;set nolabel");\r
-                   fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k2,l2,k1,l1);\r
+                   fprintf(ficgp,"\nset parametric;unset label");\r
+                   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);\r
                    fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");\r
-                   fprintf(fichtm,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup> :<a href=\"varpijgr%s%d%1d%1d-%1d%1d.png\">varpijgr%s%d%1d%1d-%1d%1d.png</A>, ",k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1);\r
-                   fprintf(fichtm,"\n<br><img src=\"varpijgr%s%d%1d%1d-%1d%1d.png\"> ",optionfilefiname, j1,k2,l2,k1,l1);\r
-                   fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\"",optionfilefiname, j1,k2,l2,k1,l1);\r
-                   fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);\r
-                   fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);\r
-                   /*              fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) t \"%d\"",\\r
-                           mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \\r
-                           mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);\r
-                   */\r
-                   fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\\r
-                           mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \\r
-                           mu1,std,v11,sqrt(lc1),v12,sqrt(lc2));\r
+                   fprintf(fichtm,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup> :<a href=\"varpijgr%s%d%1d%1d-%1d%1d.png\">varpijgr%s%d%1d%1d-%1d%1d.png</A>, ",k1,l1,k2,l2,optionfilefiname, j1,k1,l1,k2,l2,optionfilefiname, j1,k1,l1,k2,l2);\r
+                   fprintf(fichtm,"\n<br><img src=\"varpijgr%s%d%1d%1d-%1d%1d.png\"> ",optionfilefiname, j1,k1,l1,k2,l2);\r
+                   fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\"",optionfilefiname, j1,k1,l1,k2,l2);\r
+                   fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);\r
+                   fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);\r
+                   fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\\r
+                           mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),\\r
+                           mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));\r
                  }else{\r
                    first=0;\r
-                   fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);\r
-                   fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);\r
-                   /*\r
-                   fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) t \"%d\"",\\r
-                           mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \\r
-                           mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);\r
-                   */\r
-                   fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\\r
-                           mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \\r
-                           mu1,std,v11,sqrt(lc1),v12,sqrt(lc2));\r
+                   fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);\r
+                   fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);\r
+                   fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\\r
+                           mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),\\r
+                           mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));\r
                  }/* if first */\r
                } /* age mod 5 */\r
              } /* end loop age */\r
-             fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\";replot;",optionfilefiname, j1,k2,l2,k1,l1);\r
+             fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\";replot;",optionfilefiname, j1,k1,l1,k2,l2);\r
              first=1;\r
            } /*l12 */\r
          } /* k12 */\r