From 57799e1cca253eda663444432624fc44e30b5f61 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 15 Apr 2013 17:08:05 +0900 Subject: [PATCH] Adapt to new matrix transpose(), invert() and abs() methods. --- NEWS | 3 +++ examples/python/matrix_howto.py | 14 +++++++------- inst/cta/src/GCTAPointing.cpp | 2 +- inst/cta/src/GCTAResponse.cpp | 12 ++++-------- src/obs/GResponse.cpp | 6 ++---- src/opt/GOptimizerLM.cpp | 9 +++++---- src/sky/GSkyDir.cpp | 3 +-- test/test_GMatrix.cpp | 14 +++----------- test/test_GMatrixSparse.cpp | 52 +++++++++++++++++++--------------------------------- test/test_GMatrixSymmetric.cpp | 56 ++++++++++++++++++++++++-------------------------------- 10 files changed, 69 insertions(+), 102 deletions(-) diff --git a/NEWS b/NEWS index 25d709d..ae2f607 100644 --- a/NEWS +++ b/NEWS @@ -37,6 +37,9 @@ GMatrixSymmetric (formerly GSymMatrix) to column(), renamed add_col() to add_to_column(), reorder arguments of GMatrixSparse::column() and GMatrixSparse::add_to_column() methods. + The transpose(), inverse() and abs() inplace transformation methods have + been replaced by methods that return a matrix. The transpose(), + inverse(), abs(), and negate() friend functions have been removed. GModelPar The "real_" prefix has been stripped from the methods operating on diff --git a/examples/python/matrix_howto.py b/examples/python/matrix_howto.py index be0fad6..52efb0f 100755 --- a/examples/python/matrix_howto.py +++ b/examples/python/matrix_howto.py @@ -68,10 +68,10 @@ def handle_matrix(): matrix[row,column] = row + column*matrix.rows() + 1.5 # Get matrix attributes - matrix_fill = matrix.fill() # Will be 1 - matrix_minimum = matrix.min() # Will be 1.5 - matrix_maximum = matrix.max() # Will be 20.5 - matrix_sum = matrix.sum() # Will be 220.0 + matrix_fill = matrix.fill() # Will be 1 + matrix_minimum = matrix.min() # Will be 1.5 + matrix_maximum = matrix.max() # Will be 20.5 + matrix_sum = matrix.sum() # Will be 220.0 # Dump some information print("Attributes of a filled GMatrix:") @@ -89,13 +89,13 @@ def handle_matrix(): other_matrix *= 2.0 # Multiply elements by 2 other_matrix /= 2.0 # Divide elements by 2 neg_matrix = -matrix # Negate matrix - abs_matrix = gammalib.abs(neg_matrix) # Absolute value of elements + abs_matrix = neg_matrix.abs() # Absolute value of elements # And now some more tricky operations - transposed = gammalib.transpose(matrix) # First transpose matrix + transposed = matrix.transpose() # First transpose matrix mult_matrix = matrix * transposed # Now we can multiply it other_matrix *= transposed # Multiply on matrix - #inverted = gammalib.invert(matrix) # NOT YET IMPLEMENTED + #inverted = matrix.invert() # NOT YET IMPLEMENTED vector = gammalib.GVector(5) # Allocate vector multiplied = matrix * vector # Vector multiplication diff --git a/inst/cta/src/GCTAPointing.cpp b/inst/cta/src/GCTAPointing.cpp index 2a960c6..c61b7ab 100644 --- a/inst/cta/src/GCTAPointing.cpp +++ b/inst/cta/src/GCTAPointing.cpp @@ -303,7 +303,7 @@ void GCTAPointing::update(void) const Rz.eulerz(-m_dir.ra_deg()); // Compute rotation matrix - m_Rback = transpose(Ry * Rz); + m_Rback = (Ry * Rz).transpose(); // Signal that we have a valid transformation cache m_has_cache = true; diff --git a/inst/cta/src/GCTAResponse.cpp b/inst/cta/src/GCTAResponse.cpp index 0f39257..3b8f04c 100644 --- a/inst/cta/src/GCTAResponse.cpp +++ b/inst/cta/src/GCTAResponse.cpp @@ -1360,10 +1360,9 @@ double GCTAResponse::irf_diffuse(const GEvent& event, // celestial coordinates GMatrix ry; GMatrix rz; - GMatrix rot; ry.eulery(dir->dec_deg() - 90.0); rz.eulerz(-dir->ra_deg()); - rot = transpose(ry * rz); + GMatrix rot = (ry * rz).transpose(); // Setup integration kernel cta_irf_diffuse_kern_theta integrand(*this, @@ -1547,10 +1546,9 @@ double GCTAResponse::npred_radial(const GSource& source, // given by (rho,omega), into celestial coordinates. GMatrix ry; GMatrix rz; - GMatrix rot; ry.eulery(model->dec() - 90.0); rz.eulerz(-model->ra()); - rot = transpose(ry * rz); + GMatrix rot = (ry * rz).transpose(); // Compute position angle of ROI centre with respect to model // centre (radians) @@ -1726,10 +1724,9 @@ double GCTAResponse::npred_elliptical(const GSource& source, // given by (rho,omega), into celestial coordinates. GMatrix ry; GMatrix rz; - GMatrix rot; ry.eulery(model->dec() - 90.0); rz.eulerz(-model->ra()); - rot = transpose(ry * rz); + GMatrix rot = (ry * rz).transpose(); // Compute position angle of ROI centre with respect to model // centre (radians) @@ -1920,10 +1917,9 @@ double GCTAResponse::npred_diffuse(const GSource& source, // given by (theta,phi), into celestial coordinates. GMatrix ry; GMatrix rz; - GMatrix rot; ry.eulery(events->roi().centre().dec_deg() - 90.0); rz.eulerz(-events->roi().centre().ra_deg()); - rot = transpose(ry * rz); + GMatrix rot = (ry * rz).transpose(); // Setup integration kernel cta_npred_diffuse_kern_theta integrand(*this, diff --git a/src/obs/GResponse.cpp b/src/obs/GResponse.cpp index 8e73e0d..7e19628 100644 --- a/src/obs/GResponse.cpp +++ b/src/obs/GResponse.cpp @@ -418,10 +418,9 @@ double GResponse::npred_radial(const GSource& source, // by (theta,phi) into celestial coordinates. GMatrix ry; GMatrix rz; - GMatrix rot; ry.eulery(spatial->dec() - 90.0); rz.eulerz(-spatial->ra()); - rot = transpose(ry * rz); + GMatrix rot = (ry * rz).transpose(); // Set offset angle integration range double theta_min = 0.0; @@ -494,10 +493,9 @@ double GResponse::npred_elliptical(const GSource& source, // by (theta,phi) into celestial coordinates. GMatrix ry; GMatrix rz; - GMatrix rot; ry.eulery(spatial->dec() - 90.0); rz.eulerz(-spatial->ra()); - rot = transpose(ry * rz); + GMatrix rot = (ry * rz).transpose(); // Set offset angle integration range double theta_min = 0.0; diff --git a/src/opt/GOptimizerLM.cpp b/src/opt/GOptimizerLM.cpp index c7fe192..dc7fa86 100644 --- a/src/opt/GOptimizerLM.cpp +++ b/src/opt/GOptimizerLM.cpp @@ -630,8 +630,9 @@ void GOptimizerLM::iteration(GOptimizerFunction& fct, GOptimizerPars& pars) // Solve: covar * X = grad. Handle matrix problems try { - covar->cholesky_decompose(true); - *grad = covar->cholesky_solver(*grad); + *grad = covar->solve(*grad); + //GMatrixSparse decomposition = covar->cholesky_decompose(true); + //*grad = decomposition.cholesky_solver(*grad); } catch (GException::matrix_zero &e) { m_status = G_LM_SINGULAR; @@ -862,11 +863,11 @@ void GOptimizerLM::errors(GOptimizerFunction& fct, GOptimizerPars& pars) // Solve: covar * X = unit try { - covar->cholesky_decompose(true); + GMatrixSparse decomposition = covar->cholesky_decompose(true); GVector unit(npars); for (int ipar = 0; ipar < npars; ++ipar) { unit[ipar] = 1.0; - GVector x = covar->cholesky_solver(unit, true); + GVector x = decomposition.cholesky_solver(unit, true); if (x[ipar] >= 0.0) { pars.par(ipar).factor_error(sqrt(x[ipar])); } diff --git a/src/sky/GSkyDir.cpp b/src/sky/GSkyDir.cpp index 823f025..9e778bb 100644 --- a/src/sky/GSkyDir.cpp +++ b/src/sky/GSkyDir.cpp @@ -327,13 +327,12 @@ void GSkyDir::rotate_deg(const double& phi, const double& theta) // Allocate Euler and rotation matrices GMatrix ry; GMatrix rz; - GMatrix rot; // Set up rotation matrix to rotate from native coordinates to // celestial coordinates ry.eulery(m_dec * gammalib::rad2deg - 90.0); rz.eulerz(-m_ra * gammalib::rad2deg); - rot = transpose(ry * rz); + GMatrix rot = (ry * rz).transpose(); // Set up native coordinate vector double phi_rad = phi * gammalib::deg2rad; diff --git a/test/test_GMatrix.cpp b/test/test_GMatrix.cpp index e13a668..103c04c 100644 --- a/test/test_GMatrix.cpp +++ b/test/test_GMatrix.cpp @@ -393,7 +393,7 @@ void TestGMatrix::matrix_operations(void) } // Test matrix multiplication - GMatrix test3 = m_test * transpose(m_test); + GMatrix test3 = m_test * m_test.transpose(); // Check if the result matrix is as expected GMatrix ref3 = test3; @@ -588,20 +588,12 @@ void TestGMatrix::matrix_functions(void) test_value(sum, value, 1.0e-20, "Test sum function"); // Transpose function - GMatrix test1 = transpose(m_test); + GMatrix test1 = m_test.transpose(); test_assert(check_matrix(m_test), "Test source matrix"); test_assert(check_matrix_trans(test1, 1.0, 0.0), "Test transpose(GMatrix) function", "Unexpected transposed matrix:\n"+test1.print()); - // Transpose method - test1 = m_test; - test1.transpose(); - test_assert(check_matrix(m_test), "Test source matrix"); - test_assert(check_matrix_trans(test1, 1.0, 0.0), - "Test GMatrix.transpose() method", - "Unexpected transposed matrix:\n"+test1.print()); - // Convert to general matrix GMatrix test2 = GMatrix(m_test); test_assert(check_matrix(m_test), "Test source matrix"); @@ -636,7 +628,7 @@ void TestGMatrix::matrix_functions(void) } // Extract lower triangle for square matrix - GMatrix test3 = m_test * transpose(m_test); + GMatrix test3 = m_test * m_test.transpose(); GMatrix test4 = test3.extract_lower_triangle(); test_assert(check_matrix_lt(test4, test3), "Test GMatrix.extract_lower_triangle() method", diff --git a/test/test_GMatrixSparse.cpp b/test/test_GMatrixSparse.cpp index e97b093..519bbf1 100644 --- a/test/test_GMatrixSparse.cpp +++ b/test/test_GMatrixSparse.cpp @@ -708,7 +708,7 @@ void TestGMatrixSparse::matrix_operations(void) } // Test matrix multiplication - GMatrixSparse test3 = m_test * transpose(m_test); + GMatrixSparse test3 = m_test * m_test.transpose(); // Check if the result matrix is as expected GMatrixSparse ref3(g_rows, g_rows); @@ -904,21 +904,13 @@ void TestGMatrixSparse::matrix_functions(void) } test_value(sum, value, 1.0e-20, "Test sum function"); - // Transpose function - GMatrixSparse test1 = transpose(m_test); + // Transpose method + GMatrixSparse test1 = m_test.transpose(); test_assert(check_matrix(m_test), "Test source matrix"); test_assert(check_matrix_trans(test1, 1.0, 0.0), "Test transpose(GMatrixSparse) function", "Unexpected transposed matrix:\n"+test1.print()); - // Transpose method - test1 = m_test; - test1.transpose(); - test_assert(check_matrix(m_test), "Test source matrix"); - test_assert(check_matrix_trans(test1, 1.0, 0.0), - "Test GMatrixSparse.transpose() method", - "Unexpected transposed matrix:\n"+test1.print()); - // Convert to general matrix GMatrixSparse test2 = GMatrixSparse(m_test); test_assert(check_matrix(m_test), "Test source matrix"); @@ -1008,11 +1000,7 @@ void TestGMatrixSparse::matrix_cholesky(void) } // Perform Cholesky decomposition - GMatrixSparse cd = cholesky_decompose(chol_test); - - // Perform inplace Cholesky decomposition - cd = chol_test; - cd.cholesky_decompose(); + GMatrixSparse cd = chol_test.cholesky_decompose(); // Test Cholesky solver (first test) GVector e0(5); @@ -1083,8 +1071,7 @@ void TestGMatrixSparse::matrix_cholesky(void) chol_test_zero(5,5) = 1.0; // Test compressed Cholesky decomposition - GMatrixSparse cd_zero = chol_test_zero; - cd_zero.cholesky_decompose(); + GMatrixSparse cd_zero = chol_test_zero.cholesky_decompose(); // Test compressed Cholesky solver (first test) e0 = GVector(6); @@ -1155,8 +1142,7 @@ void TestGMatrixSparse::matrix_cholesky(void) chol_test_zero2(5,4) = 1.0; // Test compressed Cholesky decomposition (unsymmetric case) - GMatrixSparse cd_zero2 = chol_test_zero2; - cd_zero2.cholesky_decompose(); + GMatrixSparse cd_zero2 = chol_test_zero2.cholesky_decompose(); // Test compressed Cholesky solver (unsymmetric case) e0 = GVector(5); @@ -1222,19 +1208,20 @@ void TestGMatrixSparse::matrix_cholesky(void) unit(2,2) = 1.0; unit(3,3) = 1.0; unit(4,4) = 1.0; - GMatrixSparse chol_test_inv = chol_test; - chol_test_inv.cholesky_invert(); - GMatrixSparse ci_product = chol_test * chol_test_inv; - GMatrixSparse ci_residuals = ci_product - unit; - res = (abs(ci_residuals)).max(); - test_value(res, 0.0, 1.0e-15, "Test inplace Cholesky inverter"); + GMatrixSparse chol_test_inv = chol_test.cholesky_invert(); + GMatrixSparse ci_product = chol_test * chol_test_inv; + GMatrixSparse ci_residuals = ci_product - unit; + res = (ci_residuals.abs()).max(); + test_value(res, 0.0, 1.0e-15, "Test Cholesky inverter"); // Test Cholesky inverter - chol_test_inv = cholesky_invert(chol_test); + /* + chol_test_inv = chol_test.cholesky_invert(); ci_product = chol_test * chol_test_inv; ci_residuals = ci_product - unit; - res = (abs(ci_residuals)).max(); + res = (ci_residuals.abs()).max(); test_value(res, 0.0, 1.0e-15, "Test Cholesky inverter"); + */ // Test Cholesky inverter for compressed matrix unit = GMatrixSparse(6,6); @@ -1243,11 +1230,10 @@ void TestGMatrixSparse::matrix_cholesky(void) unit(2,2) = 1.0; unit(4,4) = 1.0; unit(5,5) = 1.0; - GMatrixSparse chol_test_zero_inv = chol_test_zero; - chol_test_zero_inv.cholesky_invert(); - GMatrixSparse ciz_product = chol_test_zero * chol_test_zero_inv; - GMatrixSparse ciz_residuals = ciz_product - unit; - res = (abs(ciz_residuals)).max(); + GMatrixSparse chol_test_zero_inv = chol_test_zero.cholesky_invert(); + GMatrixSparse ciz_product = chol_test_zero * chol_test_zero_inv; + GMatrixSparse ciz_residuals = ciz_product - unit; + res = (ciz_residuals.abs()).max(); test_value(res, 0.0, 1.0e-15, "Test compressed matrix Cholesky inverter"); // Return diff --git a/test/test_GMatrixSymmetric.cpp b/test/test_GMatrixSymmetric.cpp index c087d60..f8becd6 100644 --- a/test/test_GMatrixSymmetric.cpp +++ b/test/test_GMatrixSymmetric.cpp @@ -592,7 +592,7 @@ void TestGMatrixSymmetric::matrix_functions(void) test_value(sum, value, 1.0e-20, "Test sum function"); // Transpose function - GMatrixSymmetric test1 = transpose(m_test); + GMatrixSymmetric test1 = m_test.transpose(); test_assert(check_matrix(m_test), "Test source matrix"); test_assert(check_matrix(test1, 1.0, 0.0), "Test transpose(GMatrixSymmetric) function", @@ -659,30 +659,24 @@ void TestGMatrixSymmetric::matrix_compare(void) void TestGMatrixSymmetric::matrix_cholesky(void) { // Test Cholesky decomposition - GMatrixSymmetric cd = cholesky_decompose(m_test); - GMatrix cd_lower = cd.extract_lower_triangle(); - GMatrix cd_upper = transpose(cd_lower); - GMatrix cd_product = cd_lower * cd_upper; - GMatrix cd_residuals = GMatrix(m_test) - cd_product; - double res = (abs(cd_residuals)).max(); + GMatrixSymmetric cd = m_test.cholesky_decompose(); + GMatrix cd_lower = cd.extract_lower_triangle(); + GMatrix cd_upper = cd_lower.transpose(); + GMatrix cd_product = cd_lower * cd_upper; + GMatrix cd_residuals = GMatrix(m_test) - cd_product; + double res = (cd_residuals.abs()).max(); test_value(res, 0.0, 1.0e-15, "Test cholesky_decompose() method"); // Test compressed Cholesky decomposition GMatrixSymmetric test_zero = set_matrix_zero(); - GMatrixSymmetric cd_zero = cholesky_decompose(test_zero); - GMatrix cd_zero_lower = cd_zero.extract_lower_triangle(); - GMatrix cd_zero_upper = transpose(cd_zero_lower); - GMatrix cd_zero_product = cd_zero_lower * cd_zero_upper; - GMatrix cd_zero_residuals = GMatrix(test_zero) - cd_zero_product; - res = (abs(cd_zero_residuals)).max(); + GMatrixSymmetric cd_zero = test_zero.cholesky_decompose(); + GMatrix cd_zero_lower = cd_zero.extract_lower_triangle(); + GMatrix cd_zero_upper = cd_zero_lower.transpose(); + GMatrix cd_zero_product = cd_zero_lower * cd_zero_upper; + GMatrix cd_zero_residuals = GMatrix(test_zero) - cd_zero_product; + res = (cd_zero_residuals.abs()).max(); test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_decompose() method"); - // Test Cholesky inplace decomposition - GMatrixSymmetric test = m_test; - test.cholesky_decompose(); - GMatrix cd_lower2 = test.extract_lower_triangle(); - test_assert((cd_lower2 == cd_lower), "Test inplace cholesky_decompose() method"); - // Test Cholesky solver (first test) GVector e0(g_rows); GVector a0(g_rows); @@ -730,7 +724,7 @@ void TestGMatrixSymmetric::matrix_cholesky(void) a0[2] = 0.0; a0[3] = g_matrix[6]; s0 = cd_zero.cholesky_solver(a0) - e0; - res = max(abs(s0)); + res = max(abs(s0)); test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method"); // Test compressed Cholesky solver (second test) @@ -743,7 +737,7 @@ void TestGMatrixSymmetric::matrix_cholesky(void) a0[2] = 0.0; a0[3] = g_matrix[7]; s0 = cd_zero.cholesky_solver(a0) - e0; - res = max(abs(s0)); + res = max(abs(s0)); test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method"); // Test compressed Cholesky solver (third test) @@ -756,27 +750,25 @@ void TestGMatrixSymmetric::matrix_cholesky(void) a0[2] = 0.0; a0[3] = g_matrix[8]; s0 = cd_zero.cholesky_solver(a0) - e0; - res = max(abs(s0)); + res = max(abs(s0)); test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_solver() method"); // Test Cholesky inverter GMatrixSymmetric unit(g_rows,g_cols); unit(0,0) = unit(1,1) = unit(2,2) = 1.0; - GMatrixSymmetric test_inv = m_test; - test_inv.cholesky_invert(); - GMatrix ci_product = m_test * test_inv; - GMatrix ci_residuals = ci_product - unit; - res = (abs(ci_residuals)).max(); + GMatrixSymmetric test_inv = m_test.cholesky_invert(); + GMatrix ci_product = m_test * test_inv; + GMatrix ci_residuals = ci_product - unit; + res = (ci_residuals.abs()).max(); test_value(res, 0.0, 1.0e-15, "Test cholesky_invert method"); // Test Cholesky inverter for compressed matrix unit = GMatrixSymmetric(4,4); unit(0,0) = unit(1,1) = unit(3,3) = 1.0; - GMatrixSymmetric test_zero_inv = test_zero; - test_zero_inv.cholesky_invert(); - GMatrix ciz_product = test_zero * test_zero_inv; - GMatrix ciz_residuals = ciz_product - unit; - res = (abs(ciz_residuals)).max(); + GMatrixSymmetric test_zero_inv = test_zero.cholesky_invert(); + GMatrix ciz_product = test_zero * test_zero_inv; + GMatrix ciz_residuals = ciz_product - unit; + res = (ciz_residuals.abs()).max(); test_value(res, 0.0, 1.0e-15, "Test compressed cholesky_invert method"); // Return -- libgit2 0.21.4