Eirik Hatlen / gammalib

Commit 57799e1cca253eda663444432624fc44e30b5f61

Authored by Jürgen Knödlseder
1 parent 5b16109e

Adapt to new matrix transpose(), invert() and abs() methods.

Inline Side-by-side
Showing 10 changed files with 69 additions and 102 deletions
... ... @@ -37,6 +37,9 @@ GMatrixSymmetric (formerly GSymMatrix)
37 37 to column(), renamed add_col() to add_to_column(), reorder
38 38 arguments of GMatrixSparse::column() and GMatrixSparse::add_to_column()
39 39 methods.
  40 + The transpose(), inverse() and abs() inplace transformation methods have
  41 + been replaced by methods that return a matrix. The transpose(),
  42 + inverse(), abs(), and negate() friend functions have been removed.
40 43  
41 44 GModelPar
42 45 The "real_" prefix has been stripped from the methods operating on
... ...
examples/python/matrix_howto.py
  View file @57799e1
... ... @@ -68,10 +68,10 @@ def handle_matrix():
68 68 matrix[row,column] = row + column*matrix.rows() + 1.5
69 69  
70 70 # Get matrix attributes
71   - matrix_fill = matrix.fill() # Will be 1
72   - matrix_minimum = matrix.min() # Will be 1.5
73   - matrix_maximum = matrix.max() # Will be 20.5
74   - matrix_sum = matrix.sum() # Will be 220.0
  71 + matrix_fill = matrix.fill() # Will be 1
  72 + matrix_minimum = matrix.min() # Will be 1.5
  73 + matrix_maximum = matrix.max() # Will be 20.5
  74 + matrix_sum = matrix.sum() # Will be 220.0
75 75  
76 76 # Dump some information
77 77 print("Attributes of a filled GMatrix:")
... ... @@ -89,13 +89,13 @@ def handle_matrix():
89 89 other_matrix *= 2.0 # Multiply elements by 2
90 90 other_matrix /= 2.0 # Divide elements by 2
91 91 neg_matrix = -matrix # Negate matrix
92   - abs_matrix = gammalib.abs(neg_matrix) # Absolute value of elements
  92 + abs_matrix = neg_matrix.abs() # Absolute value of elements
93 93  
94 94 # And now some more tricky operations
95   - transposed = gammalib.transpose(matrix) # First transpose matrix
  95 + transposed = matrix.transpose() # First transpose matrix
96 96 mult_matrix = matrix * transposed # Now we can multiply it
97 97 other_matrix *= transposed # Multiply on matrix
98   - #inverted = gammalib.invert(matrix) # NOT YET IMPLEMENTED
  98 + #inverted = matrix.invert() # NOT YET IMPLEMENTED
99 99 vector = gammalib.GVector(5) # Allocate vector
100 100 multiplied = matrix * vector # Vector multiplication
101 101  
... ...
inst/cta/src/GCTAPointing.cpp
  View file @57799e1
... ... @@ -303,7 +303,7 @@ void GCTAPointing::update(void) const
303 303 Rz.eulerz(-m_dir.ra_deg());
304 304  
305 305 // Compute rotation matrix
306   - m_Rback = transpose(Ry * Rz);
  306 + m_Rback = (Ry * Rz).transpose();
307 307  
308 308 // Signal that we have a valid transformation cache
309 309 m_has_cache = true;
... ...
inst/cta/src/GCTAResponse.cpp
  View file @57799e1
... ... @@ -1360,10 +1360,9 @@ double GCTAResponse::irf_diffuse(const GEvent& event,
1360 1360 // celestial coordinates
1361 1361 GMatrix ry;
1362 1362 GMatrix rz;
1363   - GMatrix rot;
1364 1363 ry.eulery(dir->dec_deg() - 90.0);
1365 1364 rz.eulerz(-dir->ra_deg());
1366   - rot = transpose(ry * rz);
  1365 + GMatrix rot = (ry * rz).transpose();
1367 1366  
1368 1367 // Setup integration kernel
1369 1368 cta_irf_diffuse_kern_theta integrand(*this,
... ... @@ -1547,10 +1546,9 @@ double GCTAResponse::npred_radial(const GSource& source,
1547 1546 // given by (rho,omega), into celestial coordinates.
1548 1547 GMatrix ry;
1549 1548 GMatrix rz;
1550   - GMatrix rot;
1551 1549 ry.eulery(model->dec() - 90.0);
1552 1550 rz.eulerz(-model->ra());
1553   - rot = transpose(ry * rz);
  1551 + GMatrix rot = (ry * rz).transpose();
1554 1552  
1555 1553 // Compute position angle of ROI centre with respect to model
1556 1554 // centre (radians)
... ... @@ -1726,10 +1724,9 @@ double GCTAResponse::npred_elliptical(const GSource& source,
1726 1724 // given by (rho,omega), into celestial coordinates.
1727 1725 GMatrix ry;
1728 1726 GMatrix rz;
1729   - GMatrix rot;
1730 1727 ry.eulery(model->dec() - 90.0);
1731 1728 rz.eulerz(-model->ra());
1732   - rot = transpose(ry * rz);
  1729 + GMatrix rot = (ry * rz).transpose();
1733 1730  
1734 1731 // Compute position angle of ROI centre with respect to model
1735 1732 // centre (radians)
... ... @@ -1920,10 +1917,9 @@ double GCTAResponse::npred_diffuse(const GSource& source,
1920 1917 // given by (theta,phi), into celestial coordinates.
1921 1918 GMatrix ry;
1922 1919 GMatrix rz;
1923   - GMatrix rot;
1924 1920 ry.eulery(events->roi().centre().dec_deg() - 90.0);
1925 1921 rz.eulerz(-events->roi().centre().ra_deg());
1926   - rot = transpose(ry * rz);
  1922 + GMatrix rot = (ry * rz).transpose();
1927 1923  
1928 1924 // Setup integration kernel
1929 1925 cta_npred_diffuse_kern_theta integrand(*this,
... ...
src/obs/GResponse.cpp
  View file @57799e1
... ... @@ -418,10 +418,9 @@ double GResponse::npred_radial(const GSource& source,
418 418 // by (theta,phi) into celestial coordinates.
419 419 GMatrix ry;
420 420 GMatrix rz;
421   - GMatrix rot;
422 421 ry.eulery(spatial->dec() - 90.0);
423 422 rz.eulerz(-spatial->ra());
424   - rot = transpose(ry * rz);
  423 + GMatrix rot = (ry * rz).transpose();
425 424  
426 425 // Set offset angle integration range
427 426 double theta_min = 0.0;
... ... @@ -494,10 +493,9 @@ double GResponse::npred_elliptical(const GSource& source,
494 493 // by (theta,phi) into celestial coordinates.
495 494 GMatrix ry;
496 495 GMatrix rz;
497   - GMatrix rot;
498 496 ry.eulery(spatial->dec() - 90.0);
499 497 rz.eulerz(-spatial->ra());
500   - rot = transpose(ry * rz);
  498 + GMatrix rot = (ry * rz).transpose();
501 499  
502 500 // Set offset angle integration range
503 501 double theta_min = 0.0;
... ...
src/opt/GOptimizerLM.cpp
  View file @57799e1
... ... @@ -630,8 +630,9 @@ void GOptimizerLM::iteration(GOptimizerFunction& fct, GOptimizerPars& pars)
630 630  
631 631 // Solve: covar * X = grad. Handle matrix problems
632 632 try {
633   - covar->cholesky_decompose(true);
634   - *grad = covar->cholesky_solver(*grad);
  633 + *grad = covar->solve(*grad);
  634 + //GMatrixSparse decomposition = covar->cholesky_decompose(true);
  635 + //*grad = decomposition.cholesky_solver(*grad);
635 636 }
636 637 catch (GException::matrix_zero &e) {
637 638 m_status = G_LM_SINGULAR;
... ... @@ -862,11 +863,11 @@ void GOptimizerLM::errors(GOptimizerFunction& fct, GOptimizerPars& pars)
862 863  
863 864 // Solve: covar * X = unit
864 865 try {
865   - covar->cholesky_decompose(true);
  866 + GMatrixSparse decomposition = covar->cholesky_decompose(true);
866 867 GVector unit(npars);
867 868 for (int ipar = 0; ipar < npars; ++ipar) {
868 869 unit[ipar] = 1.0;
869   - GVector x = covar->cholesky_solver(unit, true);
  870 + GVector x = decomposition.cholesky_solver(unit, true);
870 871 if (x[ipar] >= 0.0) {
871 872 pars.par(ipar).factor_error(sqrt(x[ipar]));
872 873 }
... ...
src/sky/GSkyDir.cpp
  View file @57799e1
... ... @@ -327,13 +327,12 @@ void GSkyDir::rotate_deg(const double&amp; phi, const double&amp; theta)
327 327 // Allocate Euler and rotation matrices
328 328 GMatrix ry;
329 329 GMatrix rz;
330   - GMatrix rot;
331 330  
332 331 // Set up rotation matrix to rotate from native coordinates to
333 332 // celestial coordinates
334 333 ry.eulery(m_dec * gammalib::rad2deg - 90.0);
335 334 rz.eulerz(-m_ra * gammalib::rad2deg);
336   - rot = transpose(ry * rz);
  335 + GMatrix rot = (ry * rz).transpose();
337 336  
338 337 // Set up native coordinate vector
339 338 double phi_rad = phi * gammalib::deg2rad;
... ...
test/test_GMatrix.cpp
  View file @57799e1
... ... @@ -393,7 +393,7 @@ void TestGMatrix::matrix_operations(void)
393 393 }
394 394  
395 395 // Test matrix multiplication
396   - GMatrix test3 = m_test * transpose(m_test);
  396 + GMatrix test3 = m_test * m_test.transpose();
397 397  
398 398 // Check if the result matrix is as expected
399 399 GMatrix ref3 = test3;
... ... @@ -588,20 +588,12 @@ void TestGMatrix::matrix_functions(void)
588 588 test_value(sum, value, 1.0e-20, "Test sum function");
589 589  
590 590 // Transpose function
591   - GMatrix test1 = transpose(m_test);
  591 + GMatrix test1 = m_test.transpose();
592 592 test_assert(check_matrix(m_test), "Test source matrix");
593 593 test_assert(check_matrix_trans(test1, 1.0, 0.0),
594 594 "Test transpose(GMatrix) function",
595 595 "Unexpected transposed matrix:\n"+test1.print());
596 596  
597   - // Transpose method
598   - test1 = m_test;
599   - test1.transpose();
600   - test_assert(check_matrix(m_test), "Test source matrix");
601   - test_assert(check_matrix_trans(test1, 1.0, 0.0),
602   - "Test GMatrix.transpose() method",
603   - "Unexpected transposed matrix:\n"+test1.print());
604   -
605 597 // Convert to general matrix
606 598 GMatrix test2 = GMatrix(m_test);
607 599 test_assert(check_matrix(m_test), "Test source matrix");
... ... @@ -636,7 +628,7 @@ void TestGMatrix::matrix_functions(void)
636 628 }
637 629  
638 630 // Extract lower triangle for square matrix
639   - GMatrix test3 = m_test * transpose(m_test);
  631 + GMatrix test3 = m_test * m_test.transpose();
640 632 GMatrix test4 = test3.extract_lower_triangle();
641 633 test_assert(check_matrix_lt(test4, test3),
642 634 "Test GMatrix.extract_lower_triangle() method",
... ...
test/test_GMatrixSparse.cpp
  View file @57799e1
... ... @@ -708,7 +708,7 @@ void TestGMatrixSparse::matrix_operations(void)
708 708 }
709 709  
710 710 // Test matrix multiplication
711   - GMatrixSparse test3 = m_test * transpose(m_test);
  711 + GMatrixSparse test3 = m_test * m_test.transpose();
712 712  
713 713 // Check if the result matrix is as expected
714 714 GMatrixSparse ref3(g_rows, g_rows);
... ... @@ -904,21 +904,13 @@ void TestGMatrixSparse::matrix_functions(void)
904 904 }
905 905 test_value(sum, value, 1.0e-20, "Test sum function");
906 906  
907   - // Transpose function
908   - GMatrixSparse test1 = transpose(m_test);
  907 + // Transpose method
  908 + GMatrixSparse test1 = m_test.transpose();
909 909 test_assert(check_matrix(m_test), "Test source matrix");
910 910 test_assert(check_matrix_trans(test1, 1.0, 0.0),
911 911 "Test transpose(GMatrixSparse) function",
912 912 "Unexpected transposed matrix:\n"+test1.print());
913 913  
914   - // Transpose method
915   - test1 = m_test;
916   - test1.transpose();
917   - test_assert(check_matrix(m_test), "Test source matrix");
918   - test_assert(check_matrix_trans(test1, 1.0, 0.0),
919   - "Test GMatrixSparse.transpose() method",
920   - "Unexpected transposed matrix:\n"+test1.print());
921   -
922 914 // Convert to general matrix
923 915 GMatrixSparse test2 = GMatrixSparse(m_test);
924 916 test_assert(check_matrix(m_test), "Test source matrix");
... ... @@ -1008,11 +1000,7 @@ void TestGMatrixSparse::matrix_cholesky(void)
1008 1000 }
1009 1001  
1010 1002 // Perform Cholesky decomposition
1011   - GMatrixSparse cd = cholesky_decompose(chol_test);
1012   -
1013   - // Perform inplace Cholesky decomposition
1014   - cd = chol_test;
1015   - cd.cholesky_decompose();
  1003 + GMatrixSparse cd = chol_test.cholesky_decompose();
1016 1004  
1017 1005 // Test Cholesky solver (first test)
1018 1006 GVector e0(5);
... ... @@ -1083,8 +1071,7 @@ void TestGMatrixSparse::matrix_cholesky(void)
1083 1071 chol_test_zero(5,5) = 1.0;
1084 1072  
1085 1073 // Test compressed Cholesky decomposition
1086   - GMatrixSparse cd_zero = chol_test_zero;
1087   - cd_zero.cholesky_decompose();
  1074 + GMatrixSparse cd_zero = chol_test_zero.cholesky_decompose();
1088 1075  
1089 1076 // Test compressed Cholesky solver (first test)
1090 1077 e0 = GVector(6);
... ... @@ -1155,8 +1142,7 @@ void TestGMatrixSparse::matrix_cholesky(void)
1155 1142 chol_test_zero2(5,4) = 1.0;
1156 1143  
1157 1144 // Test compressed Cholesky decomposition (unsymmetric case)
1158   - GMatrixSparse cd_zero2 = chol_test_zero2;
1159   - cd_zero2.cholesky_decompose();
  1145 + GMatrixSparse cd_zero2 = chol_test_zero2.cholesky_decompose();
1160 1146  
1161 1147 // Test compressed Cholesky solver (unsymmetric case)
1162 1148 e0 = GVector(5);
... ... @@ -1222,19 +1208,20 @@ void TestGMatrixSparse::matrix_cholesky(void)
1222 1208 unit(2,2) = 1.0;
1223 1209 unit(3,3) = 1.0;
1224 1210 unit(4,4) = 1.0;
1225   - GMatrixSparse chol_test_inv = chol_test;
1226   - chol_test_inv.cholesky_invert();
1227   - GMatrixSparse ci_product = chol_test * chol_test_inv;
1228   - GMatrixSparse ci_residuals = ci_product - unit;
1229   - res = (abs(ci_residuals)).max();
1230   - test_value(res, 0.0, 1.0e-15, "Test inplace Cholesky inverter");
  1211 + GMatrixSparse chol_test_inv = chol_test.cholesky_invert();
  1212 + GMatrixSparse ci_product = chol_test * chol_test_inv;
  1213 + GMatrixSparse ci_residuals = ci_product - unit;
  1214 + res = (ci_residuals.abs()).max();
  1215 + test_value(res, 0.0, 1.0e-15, "Test Cholesky inverter");
1231 1216  
1232 1217 // Test Cholesky inverter
1233   - chol_test_inv = cholesky_invert(chol_test);
  1218 + /*
  1219 + chol_test_inv = chol_test.cholesky_invert();
1234 1220 ci_product = chol_test * chol_test_inv;
1235 1221 ci_residuals = ci_product - unit;
1236   - res = (abs(ci_residuals)).max();
  1222 + res = (ci_residuals.abs()).max();
1237 1223 test_value(res, 0.0, 1.0e-15, "Test Cholesky inverter");
  1224 + */
1238 1225  
1239 1226 // Test Cholesky inverter for compressed matrix
1240 1227 unit = GMatrixSparse(6,6);
... ... @@ -1243,11 +1230,10 @@ void TestGMatrixSparse::matrix_cholesky(void)
1243 1230 unit(2,2) = 1.0;
1244 1231 unit(4,4) = 1.0;
1245 1232 unit(5,5) = 1.0;
1246   - GMatrixSparse chol_test_zero_inv = chol_test_zero;
1247   - chol_test_zero_inv.cholesky_invert();
1248   - GMatrixSparse ciz_product = chol_test_zero * chol_test_zero_inv;
1249   - GMatrixSparse ciz_residuals = ciz_product - unit;
1250   - res = (abs(ciz_residuals)).max();
  1233 + GMatrixSparse chol_test_zero_inv = chol_test_zero.cholesky_invert();
  1234 + GMatrixSparse ciz_product = chol_test_zero * chol_test_zero_inv;
  1235 + GMatrixSparse ciz_residuals = ciz_product - unit;
  1236 + res = (ciz_residuals.abs()).max();
1251 1237 test_value(res, 0.0, 1.0e-15, "Test compressed matrix Cholesky inverter");
1252 1238  
1253 1239 // Return
... ...
test/test_GMatrixSymmetric.cpp
  View file @57799e1
... ... @@ -592,7 +592,7 @@ void TestGMatrixSymmetric::matrix_functions(void)
592 592 test_value(sum, value, 1.0e-20, "Test sum function");
593 593  
594 594 // Transpose function
595   - GMatrixSymmetric test1 = transpose(m_test);
  595 + GMatrixSymmetric test1 = m_test.transpose();
596 596 test_assert(check_matrix(m_test), "Test source matrix");
597 597 test_assert(check_matrix(test1, 1.0, 0.0),
598 598 "Test transpose(GMatrixSymmetric) function",
... ... @@ -659,30 +659,24 @@ void TestGMatrixSymmetric::matrix_compare(void)
659 659 void TestGMatrixSymmetric::matrix_cholesky(void)
660 660 {
661 661 // Test Cholesky decomposition
662   - GMatrixSymmetric cd = cholesky_decompose(m_test);
663   - GMatrix cd_lower = cd.extract_lower_triangle();
664   - GMatrix cd_upper = transpose(cd_lower);
665   - GMatrix cd_product = cd_lower * cd_upper;
666   - GMatrix cd_residuals = GMatrix(m_test) - cd_product;
667   - double res = (abs(cd_residuals)).max();
  662 + GMatrixSymmetric cd = m_test.cholesky_decompose();
  663 + GMatrix cd_lower = cd.extract_lower_triangle();
  664 + GMatrix cd_upper = cd_lower.transpose();
  665 + GMatrix cd_product = cd_lower * cd_upper;
  666 + GMatrix cd_residuals = GMatrix(m_test) - cd_product;
  667 + double res = (cd_residuals.abs()).max();
668 668 test_value(res, 0.0, 1.0e-15, "Test cholesky_decompose() method");
669 669  
670 670 // Test compressed Cholesky decomposition
671 671 GMatrixSymmetric test_zero = set_matrix_zero();
672   - GMatrixSymmetric cd_zero = cholesky_decompose(test_zero);
673   - GMatrix cd_zero_lower = cd_zero.extract_lower_triangle();
674   - GMatrix cd_zero_upper = transpose(cd_zero_lower);
675   - GMatrix cd_zero_product = cd_zero_lower * cd_zero_upper;
676   - GMatrix cd_zero_residuals = GMatrix(test_zero) - cd_zero_product;
677   - res = (abs(cd_zero_residuals)).max();
  672 + GMatrixSymmetric cd_zero = test_zero.cholesky_decompose();
  673 + GMatrix cd_zero_lower = cd_zero.extract_lower_triangle();
  674 + GMatrix cd_zero_upper = cd_zero_lower.transpose();
  675 + GMatrix cd_zero_product = cd_zero_lower * cd_zero_upper;
  676 + GMatrix cd_zero_residuals = GMatrix(test_zero) - cd_zero_product;
  677 + res = (cd_zero_residuals.abs()).max();
678 678 test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_decompose() method");
679 679  
680   - // Test Cholesky inplace decomposition
681   - GMatrixSymmetric test = m_test;
682   - test.cholesky_decompose();
683   - GMatrix cd_lower2 = test.extract_lower_triangle();
684   - test_assert((cd_lower2 == cd_lower), "Test inplace cholesky_decompose() method");
685   -
686 680 // Test Cholesky solver (first test)
687 681 GVector e0(g_rows);
688 682 GVector a0(g_rows);
... ... @@ -730,7 +724,7 @@ void TestGMatrixSymmetric::matrix_cholesky(void)
730 724 a0[2] = 0.0;
731 725 a0[3] = g_matrix[6];
732 726 s0 = cd_zero.cholesky_solver(a0) - e0;
733   - res = max(abs(s0));
  727 + res = max(abs(s0));
734 728 test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method");
735 729  
736 730 // Test compressed Cholesky solver (second test)
... ... @@ -743,7 +737,7 @@ void TestGMatrixSymmetric::matrix_cholesky(void)
743 737 a0[2] = 0.0;
744 738 a0[3] = g_matrix[7];
745 739 s0 = cd_zero.cholesky_solver(a0) - e0;
746   - res = max(abs(s0));
  740 + res = max(abs(s0));
747 741 test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method");
748 742  
749 743 // Test compressed Cholesky solver (third test)
... ... @@ -756,27 +750,25 @@ void TestGMatrixSymmetric::matrix_cholesky(void)
756 750 a0[2] = 0.0;
757 751 a0[3] = g_matrix[8];
758 752 s0 = cd_zero.cholesky_solver(a0) - e0;
759   - res = max(abs(s0));
  753 + res = max(abs(s0));
760 754 test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method");
761 755  
762 756 // Test Cholesky inverter
763 757 GMatrixSymmetric unit(g_rows,g_cols);
764 758 unit(0,0) = unit(1,1) = unit(2,2) = 1.0;
765   - GMatrixSymmetric test_inv = m_test;
766   - test_inv.cholesky_invert();
767   - GMatrix ci_product = m_test * test_inv;
768   - GMatrix ci_residuals = ci_product - unit;
769   - res = (abs(ci_residuals)).max();
  759 + GMatrixSymmetric test_inv = m_test.cholesky_invert();
  760 + GMatrix ci_product = m_test * test_inv;
  761 + GMatrix ci_residuals = ci_product - unit;
  762 + res = (ci_residuals.abs()).max();
770 763 test_value(res, 0.0, 1.0e-15, "Test cholesky_invert method");
771 764  
772 765 // Test Cholesky inverter for compressed matrix
773 766 unit = GMatrixSymmetric(4,4);
774 767 unit(0,0) = unit(1,1) = unit(3,3) = 1.0;
775   - GMatrixSymmetric test_zero_inv = test_zero;
776   - test_zero_inv.cholesky_invert();
777   - GMatrix ciz_product = test_zero * test_zero_inv;
778   - GMatrix ciz_residuals = ciz_product - unit;
779   - res = (abs(ciz_residuals)).max();
  768 + GMatrixSymmetric test_zero_inv = test_zero.cholesky_invert();
  769 + GMatrix ciz_product = test_zero * test_zero_inv;
  770 + GMatrix ciz_residuals = ciz_product - unit;
  771 + res = (ciz_residuals.abs()).max();
780 772 test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_invert method");
781 773  
782 774 // Return
... ...