diff options
author | Arnur Nigmetov <nigmetov@tugraz.at> | 2019-12-03 21:14:03 +0100 |
---|---|---|
committer | Arnur Nigmetov <nigmetov@tugraz.at> | 2020-02-18 15:02:39 +0100 |
commit | 6942d80c4d49239bca9cace9833aa74aee11ddcb (patch) | |
tree | 59797321de3b11800de28db8a4dc212828d11844 /matching/src/test_generator.cpp | |
parent | 919f04a0ad3406b9c02fbce3a9a92c03c9e85c94 (diff) |
Add matching distance code.
Diffstat (limited to 'matching/src/test_generator.cpp')
-rw-r--r-- | matching/src/test_generator.cpp | 208 |
1 files changed, 208 insertions, 0 deletions
diff --git a/matching/src/test_generator.cpp b/matching/src/test_generator.cpp new file mode 100644 index 0000000..a0d8fc7 --- /dev/null +++ b/matching/src/test_generator.cpp @@ -0,0 +1,208 @@ +#include <map> +#include <vector> +#include <random> +#include <iostream> +#include <algorithm> + +#include "opts/opts.h" +#include "spdlog/spdlog.h" +#include "spdlog/fmt/ostr.h" + +#include "common_util.h" +#include "bifiltration.h" + +using Index = md::Index; +using Point = md::Point; +using Column = md::Column; + +int g_max_coord = 100; + +using ASimplex = md::AbstractSimplex; + +using ASimplexToBirthMap = std::map<ASimplex, Point>; + +namespace spd = spdlog; + +// random generator is global +std::random_device rd; +std::mt19937_64 gen(rd()); + +//std::mt19937_64 gen(42); + +Point get_random_position(int max_coord) +{ + assert(max_coord > 0); + std::uniform_int_distribution<int> distr(0, max_coord); + return Point(distr(gen), distr(gen)); + +} + +Point get_random_position_less_than(Point ub) +{ + std::uniform_int_distribution<int> distr_x(0, ub.x); + std::uniform_int_distribution<int> distr_y(0, ub.y); + return Point(distr_x(gen), distr_y(gen)); +} + +Point get_random_position_greater_than(Point lb) +{ + std::uniform_int_distribution<int> distr_x(lb.x, g_max_coord); + std::uniform_int_distribution<int> distr_y(lb.y, g_max_coord); + return Point(distr_x(gen), distr_y(gen)); +} + +// non-proper faces (empty and simplex itself) are also faces +bool is_face(const ASimplex& face_candidate, const ASimplex& coface_candidate) +{ + return std::includes(coface_candidate.begin(), coface_candidate.end(), face_candidate.begin(), + face_candidate.end()); +} + +bool is_top_simplex(const ASimplex& candidate_simplex, const std::vector<ASimplex>& current_top_simplices) +{ + return std::none_of(current_top_simplices.begin(), current_top_simplices.end(), + [&candidate_simplex](const ASimplex& ts) { return is_face(candidate_simplex, ts); }); +} + +void add_if_top(const ASimplex& candidate_simplex, std::vector<ASimplex>& current_top_simplices) +{ + // check that candidate simplex is not face of someone in current_top_simplices + if (!is_top_simplex(candidate_simplex, current_top_simplices)) + return; + + spd::debug("candidate_simplex is top, will be added to top_simplices"); + // remove s from currrent_top_simplices, if s is face of candidate_simplex + current_top_simplices.erase(std::remove_if(current_top_simplices.begin(), current_top_simplices.end(), + [candidate_simplex](const ASimplex& s) { return is_face(s, candidate_simplex); }), + current_top_simplices.end()); + + current_top_simplices.push_back(candidate_simplex); +} + +ASimplex get_random_simplex(int n_vertices, int dim) +{ + std::vector<int> all_vertices(n_vertices, 0); + // fill in with 0..n-1 + std::iota(all_vertices.begin(), all_vertices.end(), 0); + std::shuffle(all_vertices.begin(), all_vertices.end(), gen); + return ASimplex(all_vertices.begin(), all_vertices.begin() + dim + 1, true); +} + +void generate_positions(const ASimplex& s, ASimplexToBirthMap& simplex_to_birth, Point upper_bound) +{ + auto pos = get_random_position_less_than(upper_bound); + auto curr_pos_iter = simplex_to_birth.find(s); + if (curr_pos_iter != simplex_to_birth.end()) + pos = md::greatest_lower_bound(pos, curr_pos_iter->second); + simplex_to_birth[s] = pos; + for(const ASimplex& facet : s.facets()) { + generate_positions(facet, simplex_to_birth, pos); + } +} + +md::Bifiltration get_random_bifiltration(int n_vertices, int max_dim, int n_top_simplices) +{ + ASimplexToBirthMap simplex_to_birth; + + // generate vertices + for(int i = 0; i < n_vertices; ++i) { + Point vertex_birth = get_random_position(g_max_coord / 10); + ASimplex vertex; + vertex.push_back(i); + simplex_to_birth[vertex] = vertex_birth; + } + + std::vector<ASimplex> top_simplices; + // generate top simplices + while((int)top_simplices.size() < n_top_simplices) { + std::uniform_int_distribution<int> dimension_distr(1, max_dim); + int dim = dimension_distr(gen); + auto candidate_simplex = get_random_simplex(n_vertices, dim); + spd::debug("candidate_simplex = {}", candidate_simplex); + add_if_top(candidate_simplex, top_simplices); + } + + Point upper_bound{static_cast<md::Real>(g_max_coord), static_cast<md::Real>(g_max_coord)}; + for(const auto& top_simplex : top_simplices) { + generate_positions(top_simplex, simplex_to_birth, upper_bound); + } + + std::vector<std::pair<ASimplex, Point>> simplex_birth_pairs{simplex_to_birth.begin(), simplex_to_birth.end()}; + std::vector<md::Column> boundaries{simplex_to_birth.size(), md::Column()}; + +// assign ids and save boundaries + int id = 0; + + for(int dim = 0; dim <= max_dim; ++dim) { + for(int i = 0; i < (int) simplex_birth_pairs.size(); ++i) { + ASimplex& simplex = simplex_birth_pairs[i].first; + if (simplex.dim() == dim) { + simplex.id = id++; + md::Column bdry; + for(auto& facet : simplex.facets()) { + auto facet_iter = std::find_if(simplex_birth_pairs.begin(), simplex_birth_pairs.end(), + [facet](const std::pair<ASimplex, Point>& sbp) { return facet == sbp.first; }); + assert(facet_iter != simplex_birth_pairs.end()); + assert(facet_iter->first.id >= 0); + bdry.push_back(facet_iter->first.id); + } + std::sort(bdry.begin(), bdry.end()); + boundaries[i] = bdry; + } + } + } + +// create vector of Simplex-es + std::vector<md::Simplex> simplices; + for(int i = 0; i < (int) simplex_birth_pairs.size(); ++i) { + int id = simplex_birth_pairs[i].first.id; + int dim = simplex_birth_pairs[i].first.dim(); + Point birth = simplex_birth_pairs[i].second; + Column bdry = boundaries[i]; + simplices.emplace_back(id, birth, dim, bdry); + } + +// sort by id + std::sort(simplices.begin(), simplices.end(), + [](const md::Simplex& s1, const md::Simplex& s2) { return s1.id() < s2.id(); }); + for(int i = 0; i < (int)simplices.size(); ++i) { + assert(simplices[i].id() == i); + assert(i == 0 || simplices[i].dim() >= simplices[i - 1].dim()); + } + + return md::Bifiltration(simplices.begin(), simplices.end()); +} + +int main(int argc, char** argv) +{ + spd::set_level(spd::level::info); + int n_vertices; + int max_dim; + int n_top_simplices; + std::string fname; + + using opts::Option; + using opts::PosOption; + opts::Options ops; + + bool help = false; + + ops >> Option('v', "n-vertices", n_vertices, "number of vertices") + >> Option('d', "max-dim", max_dim, "maximal dim") + >> Option('m', "max-coord", g_max_coord, "maximal coordinate") + >> Option('t', "n-top-simplices", n_top_simplices, "number of top simplices") + >> Option('h', "help", help, "show help message"); + + if (!ops.parse(argc, argv) || help || !(ops >> PosOption(fname))) { + std::cerr << "Usage: " << argv[0] << "\n" << ops << std::endl; + return 1; + } + + + auto bif1 = get_random_bifiltration(n_vertices, max_dim, n_top_simplices); + std::cout << "Generated bifiltration." << std::endl; + bif1.save(fname, md::BifiltrationFormat::rene); + std::cout << "Saved to file " << fname << std::endl; + return 0; +} + |