From 819e38d6a13e58e4f856031330e3fd58d6c1f22d Mon Sep 17 00:00:00 2001 From: Konstantinos Katrioplas Date: Thu, 17 May 2018 11:03:55 +0200 Subject: [PATCH] making sure that all matrix allocations happen at compile time when possible --- .../Optimal_bounding_box/bench_obb.cpp | 29 ++++++++++++------- .../CGAL/Optimal_bounding_box/evolution.h | 11 ++++--- .../Optimal_bounding_box/fitness_function.h | 17 ++++------- .../Optimal_bounding_box/linear_algebra.h | 26 +---------------- .../include/CGAL/Optimal_bounding_box/obb.h | 14 ++++----- .../test_optimization_algorithms.cpp | 10 +++---- .../CGAL/Eigen_linear_algebra_traits.h | 20 +++++++------ 7 files changed, 52 insertions(+), 75 deletions(-) diff --git a/Optimal_bounding_box/benchmark/Optimal_bounding_box/bench_obb.cpp b/Optimal_bounding_box/benchmark/Optimal_bounding_box/bench_obb.cpp index 85d3287be9f..fb759595d1f 100644 --- a/Optimal_bounding_box/benchmark/Optimal_bounding_box/bench_obb.cpp +++ b/Optimal_bounding_box/benchmark/Optimal_bounding_box/bench_obb.cpp @@ -1,7 +1,10 @@ #include #include +#include +#include #include #include +#include #include #include @@ -41,6 +44,7 @@ void bench_finding_obb(std::string fname) { std::ifstream input(fname); + CGAL::Eigen_linear_algebra_traits la_traits; std::vector sm_points; // import a mesh @@ -53,26 +57,29 @@ void bench_finding_obb(std::string fname) gather_mesh_points(mesh, sm_points); CGAL::Timer timer; + + // 1) measure convex hull calculation timer.start(); - - CGAL::Eigen_linear_algebra_traits la_traits; - - // 1) using convex hull - std::vector obb_points1; - CGAL::Optimal_bounding_box::find_obb(sm_points, obb_points1, la_traits, true); - + CGAL::Polyhedron_3 poly; + convex_hull_3(sm_points.begin(), sm_points.end(), poly); + std::vector ch_points(poly.points_begin(), poly.points_end()); timer.stop(); - std::cout << "found obb with convex hull: " << timer.time() << " seconds\n"; + std::cout << "takes : " << timer.time() << " seconds to find the convex hull\n"; + // 2) using convex hull timer.reset(); timer.start(); + std::vector obb_points1; + CGAL::Optimal_bounding_box::find_obb(sm_points, obb_points1, la_traits, true); + timer.stop(); + std::cout << "found obb using convex hull: " << timer.time() << " seconds\n"; - // 2) without convex hull + // 3) without convex hull + timer.reset(); + timer.start(); std::vector obb_points2; CGAL::Optimal_bounding_box::find_obb(sm_points, obb_points2, la_traits, false); - timer.stop(); - std::cout << "found obb without convex hull: " << timer.time() << " seconds\n"; timer.reset(); diff --git a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/evolution.h b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/evolution.h index 30cd901b912..dee4d14a48e 100644 --- a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/evolution.h +++ b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/evolution.h @@ -19,7 +19,6 @@ // // Author(s) : Konstantinos Katrioplas - #ifndef CGAL_EVOLUTION_H #define CGAL_EVOLUTION_H @@ -34,7 +33,7 @@ namespace Optimal_bounding_box { template class Evolution { - typedef typename Linear_algebra_traits::MatrixXd MatrixXd; + typedef typename Linear_algebra_traits::MatrixX3d MatrixXd; typedef typename Linear_algebra_traits::Matrix3d Matrix3d; typedef typename Linear_algebra_traits::Vector3d Vector3d; @@ -280,14 +279,14 @@ public: //std::cout << "pop after nelder mead: " << std::endl; //pop.show_population(); //std::cout << std::endl; - Fitness_map fitness_map_debug(population, point_data); + Fitness_map fitness_map_debug(population, point_data); Matrix3d R_now = fitness_map_debug.get_best(); std::cout << "det= " << Linear_algebra_traits::determinant(R_now) << std::endl; #endif // stopping criteria - Fitness_map fitness_map(population, point_data); - new_fit_value = fitness_map.get_best_fitness_value(point_data); + Fitness_map fitness_map(population, point_data); + new_fit_value = fitness_map.get_best_fitness_value(); double difference = new_fit_value - prev_fit_value; if(abs(difference) < tolerance * new_fit_value) stale++; @@ -301,7 +300,7 @@ public: const Matrix3d get_best() { - Fitness_map fitness_map(population, point_data); + Fitness_map fitness_map(population, point_data); return fitness_map.get_best(); } diff --git a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/fitness_function.h b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/fitness_function.h index 8f00286f0ed..981ab6bdb49 100644 --- a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/fitness_function.h +++ b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/fitness_function.h @@ -27,14 +27,11 @@ #include #include -#include - namespace CGAL { namespace Optimal_bounding_box { -template -double compute_fitness(const typename Linear_algebra_traits::Matrix3d& R, - typename Linear_algebra_traits::MatrixXd& data) +template +double compute_fitness(const Vertex& R, const Matrix& data) { // R: rotation matrix CGAL_assertion(R.cols() == 3); @@ -73,11 +70,9 @@ double compute_fitness(const typename Linear_algebra_traits::Matrix3d& R, return ((xmax - xmin) * (ymax - ymin) * (zmax - zmin)); } -template +template struct Fitness_map { - typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits; // to be added as a parameter - Fitness_map(Population& p, Matrix& points) : pop(p), points(points) {} @@ -103,14 +98,14 @@ struct Fitness_map return pop[simplex_id][vertex_id]; } - double get_best_fitness_value(Matrix& data) + double get_best_fitness_value() { const Vertex best_mat = get_best(); - return compute_fitness(best_mat, data); + return compute_fitness(best_mat, points); } Population& pop; - Matrix& points; + const Matrix& points; }; }} // end namespaces diff --git a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/linear_algebra.h b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/linear_algebra.h index 1f9dedff815..e7dbe3fc99c 100644 --- a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/linear_algebra.h +++ b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/linear_algebra.h @@ -38,26 +38,6 @@ namespace Optimal_bounding_box { typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits; #endif -/* -template -void qr_factorization(Matrix& A, Matrix& Q) -{ - Eigen::HouseholderQR qr(A); - Q = qr.householderQ(); -} - -template -void qr_factorization(std::vector& Simplex) -{ - for(int i = 0; i < Simplex.size(); ++i) - { - Eigen::HouseholderQR qr(Simplex[i]); - Matrix Q = qr.householderQ(); - Simplex[i] = Q; - } -} -*/ - template const Matrix reflection(const Matrix& S_centroid, const Matrix& S_worst) { @@ -92,19 +72,15 @@ Matrix mean(const Matrix& m1, const Matrix& m2) CGAL_assertion(m2.cols() == 3); Matrix reduction = 0.5 * m1 + 0.5 * m2; - //Matrix Q; - //reduction.qr_factorization(Q); Matrix Q = Linear_algebra_traits::qr_factorization(reduction); double det = Linear_algebra_traits::determinant(Q); return Q / det; } template -const Matrix centroid(Matrix& S1, Matrix& S2, Matrix& S3) // const? +const Matrix centroid(const Matrix& S1, const Matrix& S2, const Matrix& S3) { Matrix mean = (S1 + S2 + S3) / 3.0; - //Matrix Q; - //mean.qr_factorization(Q); Matrix Q = Linear_algebra_traits::qr_factorization(mean); double det = Linear_algebra_traits::determinant(Q); return Q / det; diff --git a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/obb.h b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/obb.h index a5d6ec79354..6b23f6350b3 100644 --- a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/obb.h +++ b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/obb.h @@ -49,8 +49,8 @@ typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits; // works on matrices only -template -void post_processing(const Matrix_dynamic& points, Vertex& R, Matrix_fixed& obb) +template +void post_processing(const Matrix& points, Vertex& R, Matrix& obb) { CGAL_assertion(points.cols() == 3); CGAL_assertion(R.rows() == 3); @@ -59,7 +59,7 @@ void post_processing(const Matrix_dynamic& points, Vertex& R, Matrix_fixed& obb) CGAL_assertion(obb.cols() == 3); // 1) rotate points with R - Matrix_dynamic rotated_points(points.rows(), points.cols()); + Matrix rotated_points(points.rows(), points.cols()); rotated_points = points * Linear_algebra_traits::transpose(R); // 2) get AABB from rotated points @@ -75,10 +75,8 @@ void post_processing(const Matrix_dynamic& points, Vertex& R, Matrix_fixed& obb) bbox = bbox_3(v_points.begin(), v_points.end()); K::Iso_cuboid_3 ic(bbox); - Matrix_fixed aabb; - // preallocate sanity: if Matrix is not preallocated in compile time - if(aabb.cols() != 3 && aabb.rows() != 8) - aabb.resize(8, 3); + // Matrix is [dynamic, 3] at compile time, so it needs allocation of rows at run time. + Matrix aabb(8, 3); for(std::size_t i = 0; i < 8; ++i) { @@ -104,7 +102,7 @@ void find_obb(std::vector& points, CGAL_assertion(obb_points.size() == 8); // eigen linear algebra traits - typedef typename LinearAlgebraTraits::MatrixXd MatrixXd; + typedef typename LinearAlgebraTraits::MatrixX3d MatrixXd; typedef typename LinearAlgebraTraits::Matrix3d Matrix3d; MatrixXd points_mat; diff --git a/Optimal_bounding_box/test/Optimal_bounding_box/test_optimization_algorithms.cpp b/Optimal_bounding_box/test/Optimal_bounding_box/test_optimization_algorithms.cpp index ae8e1879d95..763ed60fbd6 100644 --- a/Optimal_bounding_box/test/Optimal_bounding_box/test_optimization_algorithms.cpp +++ b/Optimal_bounding_box/test/Optimal_bounding_box/test_optimization_algorithms.cpp @@ -20,7 +20,7 @@ void test_nelder_mead() { typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits; typedef Linear_algebra_traits::Matrix3d Matrix3d; - typedef Linear_algebra_traits::MatrixXd MatrixXd; + typedef Linear_algebra_traits::MatrixX3d MatrixXd; MatrixXd data_points(4,3); data_points(0,0) = 0.866802; @@ -143,7 +143,7 @@ void test_nelder_mead() void test_genetic_algorithm() { - CGAL::Eigen_dense_matrix data_points(4, 3); + CGAL::Eigen_dense_matrix data_points(4, 3); // -1 : dynamic size at run time data_points(0,0) = 0.866802; data_points(0,1) = 0.740808, data_points(0,2) = 0.895304, @@ -170,7 +170,7 @@ void test_genetic_algorithm() void test_random_unit_tetra() { // this is dynamic at run times - CGAL::Eigen_dense_matrix data_points(4, 3); + CGAL::Eigen_dense_matrix data_points(4, 3); // points are on their convex hull data_points(0,0) = 0.866802; @@ -244,7 +244,7 @@ void test_reference_tetrahedron(const char* fname) } typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits; - typedef Linear_algebra_traits::MatrixXd MatrixXd; + typedef Linear_algebra_traits::MatrixX3d MatrixXd; typedef Linear_algebra_traits::Matrix3d Matrix3d; // points in a matrix @@ -279,7 +279,7 @@ void test_long_tetrahedron(std::string fname) } typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits; - typedef Linear_algebra_traits::MatrixXd MatrixXd; + typedef Linear_algebra_traits::MatrixX3d MatrixXd; typedef Linear_algebra_traits::Matrix3d Matrix3d; // points in a matrix diff --git a/Solver_interface/include/CGAL/Eigen_linear_algebra_traits.h b/Solver_interface/include/CGAL/Eigen_linear_algebra_traits.h index 5da5e294e22..c1d2716c1cd 100644 --- a/Solver_interface/include/CGAL/Eigen_linear_algebra_traits.h +++ b/Solver_interface/include/CGAL/Eigen_linear_algebra_traits.h @@ -38,7 +38,6 @@ class Eigen_dense_matrix { public: typedef Eigen::Matrix EigenType; - typedef Eigen::Matrix EigenType3; Eigen_dense_matrix(std::size_t nrows, std::size_t ncols) : m_matrix(static_cast(nrows), static_cast(ncols)) @@ -94,7 +93,6 @@ template class Eigen_dense_vector { private: - //typedef Eigen::Vector3d EigenType; typedef Eigen::Matrix EigenType; public: @@ -119,11 +117,14 @@ public: // dynamic size at run time typedef CGAL::Eigen_dense_matrix MatrixXd; + // dynamic rows in run time, fixed cols in compile time + typedef CGAL::Eigen_dense_matrix MatrixX3d; + // fixed at compile time typedef CGAL::Eigen_dense_matrix Matrix3d; // fixed at compile time - typedef CGAL::Eigen_dense_vector Vector3d; + typedef CGAL::Eigen_dense_vector Vector3d; template static Matrix transpose(const Matrix& mat) @@ -155,11 +156,11 @@ public: } template - static CGAL::Eigen_dense_vector row(CGAL::Eigen_dense_matrix& A, int i_) + static CGAL::Eigen_dense_vector row(const CGAL::Eigen_dense_matrix& A, + int i) { - return CGAL::Eigen_dense_vector(A.m_matrix.row(i_)); + return CGAL::Eigen_dense_vector(A.m_matrix.row(i)); } - }; @@ -214,7 +215,7 @@ const CGAL::Eigen_dense_matrix operator/ (const double& scalar, return CGAL::Eigen_dense_matrix (scalar / A.m_matrix); } -// addition - subtraction +// addition template const CGAL::Eigen_dense_matrix operator+ (const CGAL::Eigen_dense_matrix & A, const CGAL::Eigen_dense_matrix & B) @@ -224,9 +225,10 @@ const CGAL::Eigen_dense_matrix operator+ (const CGAL::Eigen_dense_ma // vector - matrix multiplication template -const Eigen_dense_vector operator* (const CGAL::Eigen_dense_matrix& A, CGAL::Eigen_dense_vector& V) +const Eigen_dense_vector operator* (const CGAL::Eigen_dense_matrix& A, + const CGAL::Eigen_dense_vector& V) { - return Eigen_dense_vector(A.m_matrix * V.m_vector); + return Eigen_dense_vector(A.m_matrix * V.m_vector); }