From ffedd22b13bc79e98f37f6cc76dc698b13942303 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Mon, 31 Aug 2020 12:36:30 +0200 Subject: [PATCH 01/13] initial version of BB --- include/jlpolymake/type_modules.h | 1 + src/beneath_beyond.cpp | 265 ++++++++++++++++++++++++++++++ src/jlpolymake.cpp | 2 + 3 files changed, 268 insertions(+) create mode 100644 src/beneath_beyond.cpp diff --git a/include/jlpolymake/type_modules.h b/include/jlpolymake/type_modules.h index e6e0136..0f6e94b 100644 --- a/include/jlpolymake/type_modules.h +++ b/include/jlpolymake/type_modules.h @@ -12,6 +12,7 @@ void add_array_polynomial(jlcxx::Module& jlpolymake, tparametric1 arrayt); void add_bigobject(jlcxx::Module& jlpolymake); void add_direct_calls(jlcxx::Module&); +void add_beneath_beyond(jlcxx::Module& jlpolymake); void add_incidencematrix(jlcxx::Module& jlpolymake); void add_integer(jlcxx::Module& jlpolymake); void add_matrix(jlcxx::Module& jlpolymake); diff --git a/src/beneath_beyond.cpp b/src/beneath_beyond.cpp new file mode 100644 index 0000000..5087677 --- /dev/null +++ b/src/beneath_beyond.cpp @@ -0,0 +1,265 @@ +#include "jlpolymake/jlpolymake.h" +#include +#include "jlpolymake/type_modules.h" + +namespace polymake { namespace polytope { + +template +class beneath_beyond_algo_for_ml: public beneath_beyond_algo{ + public: + typedef E value_type; + typedef const beneath_beyond_algo Base; + + beneath_beyond_algo_for_ml(): Base() + { + initialized = false; + }; + + template + void initialize(const Matrix& rays, const Matrix& lins, Iterator perm); + void initialize(const Matrix& rays, const Matrix& lins) + { + #if POLYMAKE_DEBUG + enable_debug_output(); + #endif + initialize(rays, lins, entire(sequence(0, rays.rows()))); + }; + + void process_point(Int p); + + void clear(); + + + // TODO: bundle all results in a structure, move all numbers into it + template + void compute(const Matrix& rays, const Matrix& lins, Iterator perm); + void compute(const Matrix& rays, const Matrix& lins) + { + #if POLYMAKE_DEBUG + enable_debug_output(); + #endif + compute(rays, lins, entire(sequence(0, rays.rows()))); + }; + + protected: + void stop_cleanup(); + + using Base::source_points; + using Base::source_linealities; + using Base::linealities_so_far; + using Base::expect_redundant; + using Base::source_lineality_basis; + using Base::linealities; + using Base::transform_points; + using Base::points; + using Base::generic_position; + using Base::triang_size; + using Base::AH; + using Base::interior_points; + using Base::vertices_this_step; + using Base::interior_points_this_step; + using Base::facet_normals_valid; + using Base::facet_normals_low_dim; + using Base::dual_graph; + using Base::vertices_so_far; + using Base::make_triangulation; + using Base::triangulation; + using Base::is_cone; + using Base::facets; + using compute_state = typename Base::compute_state; + using Base::state; + + class stop_calculation {}; + + private: + bool initialized; + Bitset points_added; +}; + + +template +template +void beneath_beyond_algo_for_ml::initialize(const Matrix& rays, const Matrix& lins, Iterator perm) +{ + source_points = &rays; + source_linealities = &lins; + + linealities_so_far.resize(0,rays.cols()); + + try { + if (lins.rows() != 0) { + if (expect_redundant) { + source_lineality_basis = basis_rows(lins); + linealities_so_far = lins.minor(source_lineality_basis, All); + linealities = &linealities_so_far; + } else { + linealities = source_linealities; + } + transform_points(); // the only place where stop_calculation could be thrown + } else { + points = source_points; + linealities = expect_redundant ? &linealities_so_far : source_linealities; + } + + generic_position = !expect_redundant; + triang_size = 0; + AH = unit_matrix(points->cols()); + if (expect_redundant) { + interior_points.resize(points->rows()); + vertices_this_step.resize(points->rows()); + interior_points_this_step.resize(points->rows()); + } + + state = compute_state::zero; // moved from the main compute loop + + points_added = Bitset(); + initialized = true; + } + catch (const stop_calculation&) { +#if POLYMAKE_DEBUG + if (debug >= do_dump) cout << "stop: failed to initialize beneath_beyond_algo" << endl; +#endif + // TODO: some cleanup?? + } +}; + +template +void beneath_beyond_algo_for_ml::process_point(Int p){ + if ( !points_added.contains(p) ){ + Base::process_point(p); + points_added += p; +#if POLYMAKE_DEBUG + std::cout << "processed point p = " << p << std::endl; +#endif + }; +}; + +template +template +void beneath_beyond_algo_for_ml::compute(const Matrix& rays, const Matrix& lins, Iterator perm){ + + initialize(rays, lins); + + try + { + for (; !perm.at_end(); ++perm) + process_point(*perm); + } + catch (const stop_calculation&){ +#if POLYMAKE_DEBUG + if (debug >= do_dump) cout << "stop: degenerated to full linear space" << endl; +#endif + stop_cleanup(); + } + + clear(); + +#if POLYMAKE_DEBUG + if (debug >= do_dump) { + cout << "final "; + dump(); + } +#endif + +}; + +template +void beneath_beyond_algo_for_ml::stop_cleanup(){ + state = compute_state::zero; + dual_graph.clear(); + vertices_so_far.clear(); + points = source_points; + interior_points = sequence(0, source_points->rows()); + if (make_triangulation) { + triangulation.clear(); + triang_size = 0; + } +} + +template +void beneath_beyond_algo_for_ml::clear(){ + + switch (state) { + case compute_state::zero: + if (!is_cone) { + // empty polyhedron + AH.resize(0, source_points->cols()); + linealities_so_far.resize(0, source_points->cols()); + } + break; + case compute_state::one: + // There is one empty facet in this case and the point is also a facet normal + facets[dual_graph.add_node()].normal = points->row(vertices_so_far.front()); + if (make_triangulation) { + triang_size=1; + triangulation.push_back(vertices_so_far); + } + break; + case compute_state::low_dim: + if ( !facet_normals_valid ) + { + try + { + facet_normals_low_dim(); + } + catch(const stop_calculation& ) + { + stop_cleanup(); + } + } + break; + case compute_state::full_dim: + dual_graph.squeeze(); + break; + } +} + +} +} + + +template<> struct jlcxx::IsMirroredType< + polymake::polytope::beneath_beyond_algo> : std::false_type { }; + + +void add_beneath_beyond(jlcxx::Module& polymake) +{ + polymake + .add_type>>("_BeneathBeyondAlgo") + .apply>([](auto wrapped) {}); + + polymake + .add_type>>("BeneathBeyondAlgo") + .apply>([](auto wrapped) { + typedef typename decltype(wrapped)::type WrappedT; + typedef typename decltype(wrapped)::type::value_type E; + wrapped.template constructor(); + + wrapped.method("bb_expecting_redundant", &WrappedT::expecting_redundant); + wrapped.method("bb_for_cone", &WrappedT::for_cone); + wrapped.method("bb_making_triangulation", &WrappedT::making_triangulation); + wrapped.method("bb_computing_vertices", &WrappedT::computing_vertices); + + wrapped.method("bb_compute!", static_cast< + void (polymake::polytope::beneath_beyond_algo_for_ml::*)(const pm::Matrix&, const pm::Matrix&) + >(&WrappedT::compute)); + + wrapped.method("bb_initialize!", static_cast< + void (polymake::polytope::beneath_beyond_algo_for_ml::*)(const pm::Matrix&, const pm::Matrix&) + >(&WrappedT::initialize)); + + // wrapped.method("initialize", &WrappedT::initialize); + wrapped.method("bb_add_point!", &WrappedT::process_point); + wrapped.method("bb_clear!", &WrappedT::clear); + + wrapped.method("getFacets", &WrappedT::getFacets); + wrapped.method("getVertexFacetIncidence", &WrappedT::getVertexFacetIncidence); + wrapped.method("getAffineHull", &WrappedT::getAffineHull); + wrapped.method("getVertices", &WrappedT::getVertices); + // wrapped.method("getNonRedundantPoints", &WrappedT::getNonRedundantPoints); + wrapped.method("getNonRedundantLinealities", &WrappedT::getNonRedundantLinealities); + wrapped.method("getLinealities", &WrappedT::getLinealities); + // wrapped.method("getDualGraph", &WrappedT::getDualGraph); + wrapped.method("getTriangulation", &WrappedT::getTriangulation); + }); +} diff --git a/src/jlpolymake.cpp b/src/jlpolymake.cpp index edfa0cc..fa93d6b 100644 --- a/src/jlpolymake.cpp +++ b/src/jlpolymake.cpp @@ -71,6 +71,8 @@ JLCXX_MODULE define_module_polymake(jlcxx::Module& jlpolymake) add_array_polynomial(jlpolymake, array_type); + add_beneath_beyond(jlpolymake); + add_map(jlpolymake); jlpolymake.method("initialize_polymake", &initialize_polymake); From fe9e012c388b5bda1ad6a26585f313b35eb08f04 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 11 Sep 2020 15:22:58 +0200 Subject: [PATCH 02/13] add copy constructor --- src/beneath_beyond.cpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/beneath_beyond.cpp b/src/beneath_beyond.cpp index 5087677..08ca201 100644 --- a/src/beneath_beyond.cpp +++ b/src/beneath_beyond.cpp @@ -15,6 +15,19 @@ class beneath_beyond_algo_for_ml: public beneath_beyond_algo{ initialized = false; }; + beneath_beyond_algo_for_ml(const beneath_beyond_algo_for_ml& bb) : + Base(bb), + initialized{bb.initialized}, + points_added{bb.points_added} + { + this->dual_graph.attach(this->facets); + this->dual_graph.attach(this->ridges); + if (bb.points != bb.source_points) + this->points = &(this->transformed_points); + if (bb.linealities != bb.source_linealities) + this->linealities = &(this->linealities_so_far); + }; + template void initialize(const Matrix& rays, const Matrix& lins, Iterator perm); void initialize(const Matrix& rays, const Matrix& lins) @@ -29,6 +42,10 @@ class beneath_beyond_algo_for_ml: public beneath_beyond_algo{ void clear(); + Int triangulation_size() + { + return triangulation.size(); + }; // TODO: bundle all results in a structure, move all numbers into it template From 5e4c2d6ae4016fded27aa21fad946878c2b588a0 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 11 Sep 2020 15:24:04 +0200 Subject: [PATCH 03/13] move add_beneath_beyond to jlpolymake namespace --- src/beneath_beyond.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/beneath_beyond.cpp b/src/beneath_beyond.cpp index 08ca201..499722e 100644 --- a/src/beneath_beyond.cpp +++ b/src/beneath_beyond.cpp @@ -231,13 +231,13 @@ void beneath_beyond_algo_for_ml::clear(){ } } -} -} - +} // of namespace polytope +} // of namespace polymake template<> struct jlcxx::IsMirroredType< polymake::polytope::beneath_beyond_algo> : std::false_type { }; +namespace jlpolymake { void add_beneath_beyond(jlcxx::Module& polymake) { @@ -280,3 +280,5 @@ void add_beneath_beyond(jlcxx::Module& polymake) wrapped.method("getTriangulation", &WrappedT::getTriangulation); }); } + +} // of namespace jlpolymake From 244d89f96f64dc88ea690844fbe1b0d932551d98 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 11 Sep 2020 16:42:11 +0200 Subject: [PATCH 04/13] tidy names: prefix bb_* methods with "_" --- src/beneath_beyond.cpp | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/beneath_beyond.cpp b/src/beneath_beyond.cpp index 499722e..01b8245 100644 --- a/src/beneath_beyond.cpp +++ b/src/beneath_beyond.cpp @@ -252,32 +252,32 @@ void add_beneath_beyond(jlcxx::Module& polymake) typedef typename decltype(wrapped)::type::value_type E; wrapped.template constructor(); - wrapped.method("bb_expecting_redundant", &WrappedT::expecting_redundant); - wrapped.method("bb_for_cone", &WrappedT::for_cone); - wrapped.method("bb_making_triangulation", &WrappedT::making_triangulation); - wrapped.method("bb_computing_vertices", &WrappedT::computing_vertices); + wrapped.method("_bb_expecting_redundant", &WrappedT::expecting_redundant); + wrapped.method("_bb_for_cone", &WrappedT::for_cone); + wrapped.method("_bb_making_triangulation", &WrappedT::making_triangulation); + wrapped.method("_bb_computing_vertices", &WrappedT::computing_vertices); - wrapped.method("bb_compute!", static_cast< + wrapped.method("_bb_compute!", static_cast< void (polymake::polytope::beneath_beyond_algo_for_ml::*)(const pm::Matrix&, const pm::Matrix&) >(&WrappedT::compute)); - wrapped.method("bb_initialize!", static_cast< + wrapped.method("_bb_initialize!", static_cast< void (polymake::polytope::beneath_beyond_algo_for_ml::*)(const pm::Matrix&, const pm::Matrix&) >(&WrappedT::initialize)); - // wrapped.method("initialize", &WrappedT::initialize); - wrapped.method("bb_add_point!", &WrappedT::process_point); - wrapped.method("bb_clear!", &WrappedT::clear); + wrapped.method("_bb_add_point!", &WrappedT::process_point); + wrapped.method("_bb_clear!", &WrappedT::clear); - wrapped.method("getFacets", &WrappedT::getFacets); - wrapped.method("getVertexFacetIncidence", &WrappedT::getVertexFacetIncidence); - wrapped.method("getAffineHull", &WrappedT::getAffineHull); - wrapped.method("getVertices", &WrappedT::getVertices); + wrapped.method("facets", &WrappedT::getFacets); + wrapped.method("vertex_facet_incidence", &WrappedT::getVertexFacetIncidence); + wrapped.method("affine_hull", &WrappedT::getAffineHull); + wrapped.method("vertices", &WrappedT::getVertices); // wrapped.method("getNonRedundantPoints", &WrappedT::getNonRedundantPoints); - wrapped.method("getNonRedundantLinealities", &WrappedT::getNonRedundantLinealities); - wrapped.method("getLinealities", &WrappedT::getLinealities); + wrapped.method("non_redundant_linealities", &WrappedT::getNonRedundantLinealities); + wrapped.method("linealities", &WrappedT::getLinealities); // wrapped.method("getDualGraph", &WrappedT::getDualGraph); - wrapped.method("getTriangulation", &WrappedT::getTriangulation); + wrapped.method("triangulation", &WrappedT::getTriangulation); + wrapped.method("triangulation_size", &WrappedT::triangulation_size); }); } From abe23071c3520b7fdf81dabe6d81d9217dce791f Mon Sep 17 00:00:00 2001 From: kalmarek Date: Tue, 22 Sep 2020 14:27:11 +0200 Subject: [PATCH 05/13] first try on the complete copy constructors --- src/beneath_beyond.cpp | 36 +++++++++++++++++++++++++++++++++--- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/src/beneath_beyond.cpp b/src/beneath_beyond.cpp index 01b8245..0fdc94f 100644 --- a/src/beneath_beyond.cpp +++ b/src/beneath_beyond.cpp @@ -20,12 +20,36 @@ class beneath_beyond_algo_for_ml: public beneath_beyond_algo{ initialized{bb.initialized}, points_added{bb.points_added} { + pm::cout << "copying bb" << pm::endl; this->dual_graph.attach(this->facets); this->dual_graph.attach(this->ridges); if (bb.points != bb.source_points) this->points = &(this->transformed_points); if (bb.linealities != bb.source_linealities) this->linealities = &(this->linealities_so_far); + + if (this->make_triangulation) + { + using T = typename Base::Triangulation::value_type; + std::unordered_map triangulation_map; + auto new_triangulation_iter = this->triangulation.begin(); + for (const auto& simplex : bb.triangulation) + { + triangulation_map[&simplex] = &*new_triangulation_iter; + new_triangulation_iter++; + } + + for (auto& fct_info : this->facets) + { + pm::cout << " b " << "fct_info" << pm::endl; + for (auto& s : fct_info.simplices) + { + auto new_simplex = triangulation_map[s.simplex]; + pm::cout << " replacing " << s.simplex << "with " << new_simplex << pm::endl; + s.simplex = new_simplex; + } + } + } }; template @@ -47,6 +71,11 @@ class beneath_beyond_algo_for_ml: public beneath_beyond_algo{ return triangulation.size(); }; + std::vector added_points() + { + return std::vector(points_added.begin(), points_added.end()); + } + // TODO: bundle all results in a structure, move all numbers into it template void compute(const Matrix& rays, const Matrix& lins, Iterator perm); @@ -85,10 +114,10 @@ class beneath_beyond_algo_for_ml: public beneath_beyond_algo{ using Base::facets; using compute_state = typename Base::compute_state; using Base::state; + using stop_calculation = typename Base::stop_calculation; + using Base::incident_simplex; - class stop_calculation {}; - - private: + private : bool initialized; Bitset points_added; }; @@ -278,6 +307,7 @@ void add_beneath_beyond(jlcxx::Module& polymake) // wrapped.method("getDualGraph", &WrappedT::getDualGraph); wrapped.method("triangulation", &WrappedT::getTriangulation); wrapped.method("triangulation_size", &WrappedT::triangulation_size); + wrapped.method("state", &WrappedT::added_points); }); } From 74bfc30b232393562f7350e0b6e00646e99f1af1 Mon Sep 17 00:00:00 2001 From: Benjamin Lorenz Date: Wed, 23 Sep 2020 12:42:30 +0200 Subject: [PATCH 06/13] ml_beneath_beyond: copy graph via copy_permuted to avoid creating just an alias also properly fill the NodeMap and EdgeMap --- src/beneath_beyond.cpp | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/src/beneath_beyond.cpp b/src/beneath_beyond.cpp index 0fdc94f..d217fb4 100644 --- a/src/beneath_beyond.cpp +++ b/src/beneath_beyond.cpp @@ -20,9 +20,16 @@ class beneath_beyond_algo_for_ml: public beneath_beyond_algo{ initialized{bb.initialized}, points_added{bb.points_added} { - pm::cout << "copying bb" << pm::endl; - this->dual_graph.attach(this->facets); - this->dual_graph.attach(this->ridges); + //pm::cout << "copying bb" << pm::endl; + + // we need this slightly weird copy to make sure to get a proper new graph + // instead of just an alias + Int dim = bb.dual_graph.dim(); + this->dual_graph = bb.dual_graph.copy_permuted(sequence(0,dim),sequence(0,dim)); + + this->dual_graph.attach(this->facets,entire(bb.facets)); + this->dual_graph.attach(this->ridges,entire(bb.ridges)); + if (bb.points != bb.source_points) this->points = &(this->transformed_points); if (bb.linealities != bb.source_linealities) @@ -30,26 +37,29 @@ class beneath_beyond_algo_for_ml: public beneath_beyond_algo{ if (this->make_triangulation) { + //pm::cout << "copying bb step2" << pm::endl; using T = typename Base::Triangulation::value_type; std::unordered_map triangulation_map; auto new_triangulation_iter = this->triangulation.begin(); for (const auto& simplex : bb.triangulation) { + + //pm::cout << "copying bb map: " << (void*) &simplex << " -> " << (void*) &*new_triangulation_iter << pm::endl; triangulation_map[&simplex] = &*new_triangulation_iter; new_triangulation_iter++; } - for (auto& fct_info : this->facets) { - pm::cout << " b " << "fct_info" << pm::endl; + //pm::cout << " b " << "fct_info" << pm::endl; for (auto& s : fct_info.simplices) { auto new_simplex = triangulation_map[s.simplex]; - pm::cout << " replacing " << s.simplex << "with " << new_simplex << pm::endl; + //pm::cout << " replacing " << s.simplex << "with " << new_simplex << pm::endl; s.simplex = new_simplex; } } } + //pm::cout << "copying bb done" << pm::endl; }; template @@ -116,6 +126,12 @@ class beneath_beyond_algo_for_ml: public beneath_beyond_algo{ using Base::state; using stop_calculation = typename Base::stop_calculation; using Base::incident_simplex; +#if POLYMAKE_DEBUG + using Base::enable_debug_output; + using Base::do_dump; + using Base::debug; + using Base::dump; +#endif private : bool initialized; From b94365514225221a2c3ff42453911b2204e8030e Mon Sep 17 00:00:00 2001 From: Benjamin Lorenz Date: Wed, 23 Sep 2020 13:01:10 +0200 Subject: [PATCH 07/13] ml_beneath_beyond: try making clang on macos happy --- src/beneath_beyond.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/beneath_beyond.cpp b/src/beneath_beyond.cpp index d217fb4..acb1bb6 100644 --- a/src/beneath_beyond.cpp +++ b/src/beneath_beyond.cpp @@ -125,7 +125,7 @@ class beneath_beyond_algo_for_ml: public beneath_beyond_algo{ using compute_state = typename Base::compute_state; using Base::state; using stop_calculation = typename Base::stop_calculation; - using Base::incident_simplex; + using typename Base::incident_simplex; #if POLYMAKE_DEBUG using Base::enable_debug_output; using Base::do_dump; From 1ddbb4e484b9b6d1d4806cfe2d4764b066537d4e Mon Sep 17 00:00:00 2001 From: kalmarek Date: Thu, 5 Nov 2020 16:24:09 +0100 Subject: [PATCH 08/13] =?UTF-8?q?rename=20beneath=5Fbeyond=5Falgo=5Ffor=5F?= =?UTF-8?q?ml=20=E2=86=92=20beneath=5Fbeyond=5Falgo=5Fstepwise?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/beneath_beyond.cpp | 30 +++++++++++------------------- 1 file changed, 11 insertions(+), 19 deletions(-) diff --git a/src/beneath_beyond.cpp b/src/beneath_beyond.cpp index acb1bb6..00c117a 100644 --- a/src/beneath_beyond.cpp +++ b/src/beneath_beyond.cpp @@ -5,23 +5,21 @@ namespace polymake { namespace polytope { template -class beneath_beyond_algo_for_ml: public beneath_beyond_algo{ +class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ public: typedef E value_type; typedef const beneath_beyond_algo Base; - beneath_beyond_algo_for_ml(): Base() + beneath_beyond_algo_stepwise(): Base() { initialized = false; }; - beneath_beyond_algo_for_ml(const beneath_beyond_algo_for_ml& bb) : + beneath_beyond_algo_stepwise(const beneath_beyond_algo_stepwise& bb) : Base(bb), initialized{bb.initialized}, points_added{bb.points_added} { - //pm::cout << "copying bb" << pm::endl; - // we need this slightly weird copy to make sure to get a proper new graph // instead of just an alias Int dim = bb.dual_graph.dim(); @@ -37,29 +35,23 @@ class beneath_beyond_algo_for_ml: public beneath_beyond_algo{ if (this->make_triangulation) { - //pm::cout << "copying bb step2" << pm::endl; using T = typename Base::Triangulation::value_type; std::unordered_map triangulation_map; auto new_triangulation_iter = this->triangulation.begin(); for (const auto& simplex : bb.triangulation) { - - //pm::cout << "copying bb map: " << (void*) &simplex << " -> " << (void*) &*new_triangulation_iter << pm::endl; triangulation_map[&simplex] = &*new_triangulation_iter; new_triangulation_iter++; } for (auto& fct_info : this->facets) { - //pm::cout << " b " << "fct_info" << pm::endl; for (auto& s : fct_info.simplices) { auto new_simplex = triangulation_map[s.simplex]; - //pm::cout << " replacing " << s.simplex << "with " << new_simplex << pm::endl; s.simplex = new_simplex; } } } - //pm::cout << "copying bb done" << pm::endl; }; template @@ -141,7 +133,7 @@ class beneath_beyond_algo_for_ml: public beneath_beyond_algo{ template template -void beneath_beyond_algo_for_ml::initialize(const Matrix& rays, const Matrix& lins, Iterator perm) +void beneath_beyond_algo_stepwise::initialize(const Matrix& rays, const Matrix& lins, Iterator perm) { source_points = &rays; source_linealities = &lins; @@ -186,7 +178,7 @@ void beneath_beyond_algo_for_ml::initialize(const Matrix& rays, const Matr }; template -void beneath_beyond_algo_for_ml::process_point(Int p){ +void beneath_beyond_algo_stepwise::process_point(Int p){ if ( !points_added.contains(p) ){ Base::process_point(p); points_added += p; @@ -198,7 +190,7 @@ void beneath_beyond_algo_for_ml::process_point(Int p){ template template -void beneath_beyond_algo_for_ml::compute(const Matrix& rays, const Matrix& lins, Iterator perm){ +void beneath_beyond_algo_stepwise::compute(const Matrix& rays, const Matrix& lins, Iterator perm){ initialize(rays, lins); @@ -226,7 +218,7 @@ void beneath_beyond_algo_for_ml::compute(const Matrix& rays, const Matrix< }; template -void beneath_beyond_algo_for_ml::stop_cleanup(){ +void beneath_beyond_algo_stepwise::stop_cleanup(){ state = compute_state::zero; dual_graph.clear(); vertices_so_far.clear(); @@ -239,7 +231,7 @@ void beneath_beyond_algo_for_ml::stop_cleanup(){ } template -void beneath_beyond_algo_for_ml::clear(){ +void beneath_beyond_algo_stepwise::clear(){ switch (state) { case compute_state::zero: @@ -292,7 +284,7 @@ void add_beneath_beyond(jlcxx::Module& polymake) polymake .add_type>>("BeneathBeyondAlgo") - .apply>([](auto wrapped) { + .apply>([](auto wrapped) { typedef typename decltype(wrapped)::type WrappedT; typedef typename decltype(wrapped)::type::value_type E; wrapped.template constructor(); @@ -303,11 +295,11 @@ void add_beneath_beyond(jlcxx::Module& polymake) wrapped.method("_bb_computing_vertices", &WrappedT::computing_vertices); wrapped.method("_bb_compute!", static_cast< - void (polymake::polytope::beneath_beyond_algo_for_ml::*)(const pm::Matrix&, const pm::Matrix&) + void (polymake::polytope::beneath_beyond_algo_stepwise::*)(const pm::Matrix&, const pm::Matrix&) >(&WrappedT::compute)); wrapped.method("_bb_initialize!", static_cast< - void (polymake::polytope::beneath_beyond_algo_for_ml::*)(const pm::Matrix&, const pm::Matrix&) + void (polymake::polytope::beneath_beyond_algo_stepwise::*)(const pm::Matrix&, const pm::Matrix&) >(&WrappedT::initialize)); wrapped.method("_bb_add_point!", &WrappedT::process_point); From d3625eb9fad7a916238d3eac1ef417bc91914727 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 11 Nov 2020 15:10:53 +0100 Subject: [PATCH 09/13] =?UTF-8?q?rename=20clear=20=E2=86=92=20finish?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/beneath_beyond.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/beneath_beyond.cpp b/src/beneath_beyond.cpp index 00c117a..52c1ebe 100644 --- a/src/beneath_beyond.cpp +++ b/src/beneath_beyond.cpp @@ -66,7 +66,7 @@ class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ void process_point(Int p); - void clear(); + void finish(); Int triangulation_size() { @@ -206,7 +206,7 @@ void beneath_beyond_algo_stepwise::compute(const Matrix& rays, const Matri stop_cleanup(); } - clear(); + finish(); #if POLYMAKE_DEBUG if (debug >= do_dump) { @@ -231,7 +231,7 @@ void beneath_beyond_algo_stepwise::stop_cleanup(){ } template -void beneath_beyond_algo_stepwise::clear(){ +void beneath_beyond_algo_stepwise::finish(){ switch (state) { case compute_state::zero: @@ -303,7 +303,7 @@ void add_beneath_beyond(jlcxx::Module& polymake) >(&WrappedT::initialize)); wrapped.method("_bb_add_point!", &WrappedT::process_point); - wrapped.method("_bb_clear!", &WrappedT::clear); + wrapped.method("_bb_finish!", &WrappedT::finish); wrapped.method("facets", &WrappedT::getFacets); wrapped.method("vertex_facet_incidence", &WrappedT::getVertexFacetIncidence); From f5fbb1cfce7d296c5032dd47973f951c905bf36a Mon Sep 17 00:00:00 2001 From: kalmarek Date: Thu, 12 Nov 2020 11:11:15 +0100 Subject: [PATCH 10/13] provisional deepcopy of beneath_beyond_algo_stepwise --- src/beneath_beyond.cpp | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/src/beneath_beyond.cpp b/src/beneath_beyond.cpp index 52c1ebe..db2caaa 100644 --- a/src/beneath_beyond.cpp +++ b/src/beneath_beyond.cpp @@ -4,6 +4,16 @@ namespace polymake { namespace polytope { +template +T deepcopy(const T& x){ + return T{entire(x)}; +} + +template +Vector deepcopy(const Vector& v){ + return Vector(v.dim(), entire(v)); +} + template class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ public: @@ -17,9 +27,19 @@ class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ beneath_beyond_algo_stepwise(const beneath_beyond_algo_stepwise& bb) : Base(bb), + // TODO: + // dual_graph + // AH + // facet_nullspace + initialized{bb.initialized}, points_added{bb.points_added} { + this->interior_points = deepcopy(bb.interior_points); + this->source_lineality_basis = deepcopy(bb.source_lineality_basis); + this->points_in_lineality_basis = deepcopy(bb.points_in_lineality_basis); + this->vertices_so_far = deepcopy(bb.vertices_so_far); + // we need this slightly weird copy to make sure to get a proper new graph // instead of just an alias Int dim = bb.dual_graph.dim(); @@ -33,6 +53,11 @@ class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ if (bb.linealities != bb.source_linealities) this->linealities = &(this->linealities_so_far); + for (auto ridge = entire(this->ridges); !ridge.at_end(); ++ridge) + { + *ridge = deepcopy(*ridge); + } + if (this->make_triangulation) { using T = typename Base::Triangulation::value_type; @@ -40,11 +65,15 @@ class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ auto new_triangulation_iter = this->triangulation.begin(); for (const auto& simplex : bb.triangulation) { + auto new_simplex = deepcopy(simplex); + *new_triangulation_iter = new_simplex; triangulation_map[&simplex] = &*new_triangulation_iter; new_triangulation_iter++; } for (auto& fct_info : this->facets) { + fct_info.normal = deepcopy(fct_info.normal); + fct_info.vertices = deepcopy(fct_info.vertices); for (auto& s : fct_info.simplices) { auto new_simplex = triangulation_map[s.simplex]; From 90ea63ecd1d96ddd8da9bc3adae71065127c2244 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Mon, 16 Nov 2020 11:24:39 +0100 Subject: [PATCH 11/13] deepcopy bb.AH and bb.facet_nullspace --- src/beneath_beyond.cpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/beneath_beyond.cpp b/src/beneath_beyond.cpp index db2caaa..17a2531 100644 --- a/src/beneath_beyond.cpp +++ b/src/beneath_beyond.cpp @@ -14,6 +14,12 @@ Vector deepcopy(const Vector& v){ return Vector(v.dim(), entire(v)); } +template +ListMatrix deepcopy(const ListMatrix& m){ + return ListMatrix(Matrix(m)); +} + + template class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ public: @@ -27,14 +33,11 @@ class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ beneath_beyond_algo_stepwise(const beneath_beyond_algo_stepwise& bb) : Base(bb), - // TODO: - // dual_graph - // AH - // facet_nullspace - initialized{bb.initialized}, points_added{bb.points_added} { + this->AH = deepcopy(bb.AH); + this->facet_nullspace = deepcopy(bb.facet_nullspace); this->interior_points = deepcopy(bb.interior_points); this->source_lineality_basis = deepcopy(bb.source_lineality_basis); this->points_in_lineality_basis = deepcopy(bb.points_in_lineality_basis); From 3e9d8983c65f0209787c7ac7b65e072ba39bc4f5 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Mon, 16 Nov 2020 11:25:39 +0100 Subject: [PATCH 12/13] explicitly copy Bitsets and deque used in per-step operation --- src/beneath_beyond.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/beneath_beyond.cpp b/src/beneath_beyond.cpp index 17a2531..6b36c6e 100644 --- a/src/beneath_beyond.cpp +++ b/src/beneath_beyond.cpp @@ -41,6 +41,11 @@ class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ this->interior_points = deepcopy(bb.interior_points); this->source_lineality_basis = deepcopy(bb.source_lineality_basis); this->points_in_lineality_basis = deepcopy(bb.points_in_lineality_basis); + this->vertices_this_step = Bitset{}; + this->interior_points_this_step = Bitset{}; + this->visited_facets = Bitset{}; + this->facet_queue = std::deque{}; + this->vertices_so_far = deepcopy(bb.vertices_so_far); // we need this slightly weird copy to make sure to get a proper new graph From 5b5dd4404733027d98e2eb980e78169f27613e8f Mon Sep 17 00:00:00 2001 From: Benjamin Lorenz Date: Tue, 24 Nov 2020 15:10:12 +0100 Subject: [PATCH 13/13] more copies, store full points as well --- src/beneath_beyond.cpp | 38 ++++++++++++++++++++++++++++++-------- 1 file changed, 30 insertions(+), 8 deletions(-) diff --git a/src/beneath_beyond.cpp b/src/beneath_beyond.cpp index 6b36c6e..0c06a1d 100644 --- a/src/beneath_beyond.cpp +++ b/src/beneath_beyond.cpp @@ -19,6 +19,11 @@ ListMatrix deepcopy(const ListMatrix& m){ return ListMatrix(Matrix(m)); } +template +Matrix deepcopy(const Matrix& m) { + return Matrix(m.rows(),m.cols(),entire(concat_rows(m))); +} + template class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ @@ -34,7 +39,9 @@ class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ beneath_beyond_algo_stepwise(const beneath_beyond_algo_stepwise& bb) : Base(bb), initialized{bb.initialized}, - points_added{bb.points_added} + points_added{bb.points_added}, + full_points{deepcopy(bb.full_points)}, + full_linealities{deepcopy(bb.full_linealities)} { this->AH = deepcopy(bb.AH); this->facet_nullspace = deepcopy(bb.facet_nullspace); @@ -46,6 +53,10 @@ class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ this->visited_facets = Bitset{}; this->facet_queue = std::deque{}; + this->transformed_points = deepcopy(bb.transformed_points); + this->linealities_so_far = deepcopy(bb.linealities_so_far); + this->lineality_transform = deepcopy(bb.lineality_transform); + this->vertices_so_far = deepcopy(bb.vertices_so_far); // we need this slightly weird copy to make sure to get a proper new graph @@ -56,10 +67,17 @@ class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ this->dual_graph.attach(this->facets,entire(bb.facets)); this->dual_graph.attach(this->ridges,entire(bb.ridges)); + this->source_points = &(this->full_points); if (bb.points != bb.source_points) this->points = &(this->transformed_points); + else + this->points = &(this->full_points); + + this->source_linealities = &(this->full_linealities); if (bb.linealities != bb.source_linealities) this->linealities = &(this->linealities_so_far); + else + this->linealities = &(this->full_linealities); for (auto ridge = entire(this->ridges); !ridge.at_end(); ++ridge) { @@ -165,6 +183,8 @@ class beneath_beyond_algo_stepwise: public beneath_beyond_algo{ private : bool initialized; Bitset points_added; + Matrix full_points; + Matrix full_linealities; }; @@ -172,16 +192,18 @@ template template void beneath_beyond_algo_stepwise::initialize(const Matrix& rays, const Matrix& lins, Iterator perm) { - source_points = &rays; - source_linealities = &lins; + full_points = deepcopy(rays); + full_linealities = deepcopy(lins); + source_points = &full_points; + source_linealities = &full_linealities; - linealities_so_far.resize(0,rays.cols()); + linealities_so_far.resize(0,full_points.cols()); try { - if (lins.rows() != 0) { + if (full_linealities.rows() != 0) { if (expect_redundant) { - source_lineality_basis = basis_rows(lins); - linealities_so_far = lins.minor(source_lineality_basis, All); + source_lineality_basis = basis_rows(full_linealities); + linealities_so_far = full_linealities.minor(source_lineality_basis, All); linealities = &linealities_so_far; } else { linealities = source_linealities; @@ -220,7 +242,7 @@ void beneath_beyond_algo_stepwise::process_point(Int p){ Base::process_point(p); points_added += p; #if POLYMAKE_DEBUG - std::cout << "processed point p = " << p << std::endl; + if (debug >= do_dump) std::cout << "processed point p = " << p << std::endl; #endif }; };