From b6a080a37d58c0f7e5b9354f0545400991a19fb6 Mon Sep 17 00:00:00 2001 From: crvs Date: Mon, 6 Feb 2017 18:09:05 +0100 Subject: [PATCH 01/24] finished the necessary facilities to calculate paths --- CMakeLists.txt | 4 +- lib/scomplex/graph_utils.hpp | 66 ++++++++++++++-------------- lib/scomplex/nn_utils.hpp | 78 +++++++++++++++++++++++++--------- lib/scomplex/path_snapper.hpp | 32 ++++++++++++++ lib/scomplex/qhull_parsing.hpp | 6 ++- src/libtest.cpp | 19 ++++++--- src/nn_test.cpp | 4 ++ 7 files changed, 149 insertions(+), 60 deletions(-) create mode 100644 lib/scomplex/path_snapper.hpp create mode 100644 src/nn_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 143319a..7eeadf1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -34,6 +34,8 @@ set(YOTHER_STUPID_SOURCES "./src/graphtest.cpp") add_executable(yother_stupid_test ${YOTHER_STUPID_SOURCES}) # set(CMAKE_CXX_FLAGS "-lspatialindex_c -lspatialindex ${CMAKE_CXX_FLAGS}") - set(OTHER_STUPID_SOURCES "./src/libtest.cpp") +set(OTHER_STUPID_SOURCES "./src/libtest.cpp") add_executable(my_other_stupid_test ${OTHER_STUPID_SOURCES}) +add_executable(nn_test "./src/nn_test") + diff --git a/lib/scomplex/graph_utils.hpp b/lib/scomplex/graph_utils.hpp index 3cfd5f1..2672742 100644 --- a/lib/scomplex/graph_utils.hpp +++ b/lib/scomplex/graph_utils.hpp @@ -19,7 +19,7 @@ typedef typename boost::adjacency_list< // boost::vecS, // boost::vecS, // boost::undirectedS, // - boost::no_property, // + boost::no_property, // vertex property boost::property> // Graph; @@ -61,51 +61,55 @@ Graph calculate_one_skelleton_graph( return one_skelleton; } -// calculate shortest paths: +/** + * @brief find the shortest path between two vertices in a graph. + * + * @param g: graph (passed by ref) + * @param s: source vertex (type: decltype(g)::vertex_descriptor) + * @param t: source vertex (type: decltype(g)::vertex_descriptor) + * + * @returns: a vector of vertices + * (type: std::vector) + */ std::vector shortest_path(const Graph& g, vertex_t s, vertex_t t) { + // make a predecessor map boost::vector_property_map predecessors(boost::num_vertices(g)); - boost::dijkstra_shortest_paths(g, s, boost::predecessor_map(&predecessors[0])); - + // construct the path (going backwards from the target to the source) vertex_t it = t; std::vector s_t_path{}; - while (it != s) { - s_t_path.push_back(it); + // WARNING: the push_back and reverse strategy may be more efficient + // + // use a do-while block, we want the guard to be checked afterwards + do { + // push to the front + s_t_path.insert(s_t_path.begin(), it); it = predecessors[it]; - } - s_t_path.push_back(it); - std::reverse(s_t_path.begin(), s_t_path.end()); + } while (it != s); return s_t_path; - - // std::cout << "I'm calculating shortest paths\n"; } std::vector complete_path(const Graph& g, std::vector vec) { std::vector::iterator s{vec.begin()}; - std::list full_path{*s}; - while (s != vec.end()) { - std::vector::iterator t = std::next(s); - // calculate the path between s and t - std::vector path_portion{shortest_path(g, *s, *t)}; + std::vector::iterator t = std::next(s); - // add the path between s and t to the full path + // start calculatinng the full path + std::vector full_path; + full_path.push_back(*s); + while (t != vec.end()) { + // get the next vertex + std::cout << "looking at edge" << *s << ", " << *t << "\n"; + + // calculate and add the path between s and t (contains vertices other + // than s) to the full path + std::vector path_portion{shortest_path(g, *s, *t)}; for (vertex_t vert : path_portion) { - // allways skip the first element since we are calculating the paths - // between s,t then s',t' where s' = t and that would make use of a - // non-existing edge. - if (vert != *path_portion.begin()) { - full_path.push_back(vert); - } + full_path.push_back(vert); } - + std::cout << "\n"; std::advance(s, 1); + std::advance(t, 1); } - // initialize with size so we don't have to relocate - std::vector full_path_vect(full_path.size()); - for (vertex_t vert : full_path) { - full_path_vect.push_back(vert); - } - - return full_path_vect; + return full_path; } diff --git a/lib/scomplex/nn_utils.hpp b/lib/scomplex/nn_utils.hpp index fb42161..b702403 100644 --- a/lib/scomplex/nn_utils.hpp +++ b/lib/scomplex/nn_utils.hpp @@ -5,42 +5,78 @@ #include -// #include #include +#include #include #include -// #include #include +#include typedef std::vector point_t; -typedef CGAL::Cartesian_d cart; -typedef cart::Point_d point_d; -typedef CGAL::Search_traits_d traits; -typedef CGAL::Orthogonal_k_neighbor_search neighbor_search_t; -typedef neighbor_search_t::Tree tree_t; +// kd-tree point type +typedef CGAL::Cartesian_d::Point_d point_d; // d-dimensional point -// typedef CGAL::Kd_tree kd_tree; -// typedef CGAL::K_neighbor_search k_neighbor; -// typedef CGAL::Euclidean_distance euclidean; +typedef CGAL::Search_traits_d> cartesian_traits; -// template -// typedef cart::Point_d point_d; +typedef CGAL::Search_traits_adapter< // define the search tratis + boost::tuple, // pairing d-point with its index + CGAL::Nth_of_tuple_property_map< // custom kernel + 0, // position of the points to compare + boost::tuple>, // type tha will be used + cartesian_traits> // traits to use + traits; // name -void make_tree(tree_t& tree, std::list& point_list) { - tree.reserve(point_list.size()); - std::list insertable_points; - for (auto pt : point_list) { +// type of search tree +typedef CGAL::Orthogonal_k_neighbor_search neighbor_search_t; + +typedef neighbor_search_t::Tree tree_t; + +void make_tree(tree_t& tree, std::vector& point_list) { + // tree.reserve(point_list.size()); + for (size_t ind = 0; ind < point_list.size(); ++ind) { + auto pt = point_list.at(ind); point_d pt_d(pt.size(), pt.begin(), pt.end()); - tree.insert(pt_d); + tree.insert(boost::make_tuple(pt_d, ind)); } - // return tree; -}; +} -point_t nearest_neighbour(tree_t& tree, point_t point) { +// given a point snap it to the tree +point_t nearest_neighbor(tree_t& tree, point_t point) { point_d pt(point.size(), point.begin(), point.end()); + // "1" refers to the number of nearest neighbors to find neighbor_search_t search(tree, pt, 1); - point_d result_pt = search.begin()->first; + point_d result_pt; + size_t ind; + boost::tie(result_pt, ind) = search.begin()->first; point_t point_vec(result_pt.cartesian_begin(), result_pt.cartesian_end()); return point_vec; } + +// given a vector of points snap them to the tree +std::vector snap_points(tree_t& tree, std::vector points) { + std::vector snapped_points; + for (auto point : points) { + snapped_points.push_back(nearest_neighbor(tree, point)); + } + return snapped_points; +} + +size_t nearest_neighbor_index(tree_t& tree, point_t point) { + point_d pt(point.size(), point.begin(), point.end()); + // "1" refers to the number of nearest neighbors to find + neighbor_search_t search(tree, pt, 1); + point_d result_pt; + size_t ind; + boost::tie(result_pt, ind) = search.begin()->first; + return ind; +} + +std::vector snap_points_to_indexes(tree_t& tree, + std::vector points) { + std::vector snapped_points; + for (auto point : points) { + snapped_points.push_back(nearest_neighbor_index(tree, point)); + } + return snapped_points; +} diff --git a/lib/scomplex/path_snapper.hpp b/lib/scomplex/path_snapper.hpp new file mode 100644 index 0000000..740aa4a --- /dev/null +++ b/lib/scomplex/path_snapper.hpp @@ -0,0 +1,32 @@ +#pragma once + +namespace snap { + +#include +#include +#include + +class path_snapper { + private: + tree_t point_tree; // defined in nn_utils.hpp + Graph vertex_graph; // defined in graph_utils.hpp + + public: + path_snapper(simplicial::SimplicialComplex& sc) { + make_tree(point_tree, sc.points); + vertex_graph = calculate_one_skelleton_graph(sc); + } + + std::vector snap_path(std::vector path) { + std::vector way_points = + snap_points_to_indexes(point_tree, path); + std::copy(way_points.begin(), way_points.end(), + std::ostream_iterator(std::cout, " ")); + std::cout << '\n'; + + std::vector snapped_path = + complete_path(vertex_graph, way_points); + return snapped_path; + }; +}; +}; diff --git a/lib/scomplex/qhull_parsing.hpp b/lib/scomplex/qhull_parsing.hpp index 7a88b97..fdd384b 100644 --- a/lib/scomplex/qhull_parsing.hpp +++ b/lib/scomplex/qhull_parsing.hpp @@ -1,4 +1,3 @@ - #pragma once #include @@ -32,6 +31,11 @@ void output_vector(std::vector vec) { std::pair, std::vector> parse_qhull_file( std::string filename) { std::ifstream file_stream(filename); + + if (!file_stream.is_open()) { + throw std::runtime_error("failed to open \"" + filename + "\"\n"); + } + std::string line; int current_line_number{0}; diff --git a/src/libtest.cpp b/src/libtest.cpp index 9a7a696..e3e5b64 100644 --- a/src/libtest.cpp +++ b/src/libtest.cpp @@ -1,8 +1,9 @@ #include #include -#include +//#include #include #include +#include #include "boost/graph/graph_traits.hpp" @@ -128,18 +129,17 @@ int main() { typedef typename decltype(G)::vertex_descriptor vertex_t; - /* std::cout << '\n'; std::copy(p.begin(), p.end(), std::ostream_iterator{std::cout, " "}); std::cout << '\n'; - */ + std::vector point_vec{point_list.begin(), point_list.end()}; tree_t rt; - make_tree(rt, point_list); + make_tree(rt, point_vec); point_t q{0.1, 0.2}; - point_t q_prime = nearest_neighbour(rt, q); + point_t q_prime = nearest_neighbor(rt, q); // get_r_tree_from_points<2>(point_list); std::ostream_iterator outstr(std::cout, " "); std::copy(q_prime.begin(), q_prime.end(), outstr); @@ -147,7 +147,7 @@ int main() { std::vector points_v; std::vector cells_v; std::tie(points_v, cells_v) = - parse_qhull_file("/home/crvs/qhull-test/qh-test.dat"); + parse_qhull_file("/home/crvs/dev-cpp/examples/qhull-test/qh-test.dat"); std::list points(points_v.begin(), points_v.end()); std::list> cells; for (auto cell_v : cells_v) { @@ -156,5 +156,12 @@ int main() { } simplicial::SimplicialComplex rsc(points, cells); + + snap::path_snapper snapper(rsc); + + std::cout << "\n"; + auto path = snapper.snap_path({{-1, -1}, {1, 1}}); + std::copy(path.begin(), path.end(), + std::ostream_iterator(std::cout, " ")); return 0; } diff --git a/src/nn_test.cpp b/src/nn_test.cpp new file mode 100644 index 0000000..f6fb1f7 --- /dev/null +++ b/src/nn_test.cpp @@ -0,0 +1,4 @@ +#include +#include + +int main() { return 0; } From 7ce7bd6bcc3c2aead6967af964601f1b347853bd Mon Sep 17 00:00:00 2001 From: crvs Date: Tue, 7 Feb 2017 15:20:35 +0100 Subject: [PATCH 02/24] updating utils --- lib/scomplex/Simplicial_Complex.h | 83 +++++++------------------------ lib/scomplex/base_utils.hpp | 9 ++-- lib/scomplex/graph_utils.hpp | 3 +- lib/scomplex/nn_utils.hpp | 14 +++--- lib/scomplex/path_snapper.hpp | 24 ++++++++- 5 files changed, 54 insertions(+), 79 deletions(-) diff --git a/lib/scomplex/Simplicial_Complex.h b/lib/scomplex/Simplicial_Complex.h index 780608c..ccfab6e 100644 --- a/lib/scomplex/Simplicial_Complex.h +++ b/lib/scomplex/Simplicial_Complex.h @@ -22,39 +22,31 @@ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * */ -#ifndef QUOT_SIMPLICIAL_COMPLEX -#define QUOT_SIMPLICIAL_COMPLEX +#pragma once #include -#include #include +#include #include -// #include -#include -#include #include - -#include - -// point_t needs to have an empty constructor, and a size() function call, -// further it needs a begin() function call. Mostly it needs to behave like a -// std::vector +#include +#include namespace simplicial { -typedef std::pair> chain_t; template class SimplicialComplex { private: // the type of matrix to be used typedef typename Eigen::SparseMatrix matrix_t; + typedef typename Eigen::SparseVector chain_t; + // simplex tree options struct SimpleOptions : Gudhi::Simplex_tree_options_full_featured { - // simplex tree options - typedef int Vertex_handle; + typedef size_t Vertex_handle; }; typedef typename Gudhi::Simplex_tree SimplexTree; typedef typename Gudhi::Simplex_tree::Simplex_handle @@ -62,13 +54,14 @@ class SimplicialComplex { typedef std::vector level_t; typedef std::vector levels_t; + typedef std::vector simplex_t; public: - std::vector> get_level(int dimen) { - std::vector> level; + std::vector get_level(int dimen) { + std::vector level; for (auto simp : simplices.complex_simplex_range()) { if (simplices.dimension(simp) == dimen) { - std::vector v_simp; + simplex_t v_simp; for (auto v : simplices.simplex_vertex_range(simp)) { v_simp.push_back(v); } @@ -105,8 +98,6 @@ class SimplicialComplex { if (simplices.dimension() < d) { simplices.set_dimension(d); } - // error (something is out-of-range) - // calculate_matrices(); } // we need to count i--simplices in order to get a key for them @@ -207,58 +198,19 @@ class SimplicialComplex { // need to do the boundary inclusion crap (you know what I mean) int boundary_index(Simplex_handle s_1, Simplex_handle s_2) { - // s1 has to include into s2 - - /* - int dim1 = simplices.dimension(s_1); - int dim2 = simplices.dimension(s_2); - decltype(s_1) bdry_h; - decltype(s_2) simp_h; - if (dim2 == dim1 - 1) { - simp_h = decltype(s_1)(s_1); - bdry_h = decltype(s_2)(s_2); - } else if (dim1 == dim2 - 1) { - simp_h = decltype(s_2)(s_2); - bdry_h = decltype(s_1)(s_1); - } else { - std::cerr << "CAN'T BE A BOUNDARY RELATION"; - throw; - } - - bool found = false; - int index = dim1; - - // <--- this is not working need to find the parent!! - // simp_h = simp_h.parent(); - - auto someth = simp_h->second.children()->parent(); - auto someth2 = someth.first; - - std::cout << "simp_h: " << typeid(simp_h).name() << std::endl; - std::cout << "simp_h-fst: " << someth2 << std::endl; - std::cout << "someth: " << typeid(someth).name() << std::endl; - // std::cout << "someth-fst: " << someth2 << std::endl; - - while (not found) { // <-- this is not working yet - if (simp_h == bdry_h) {found = true; index--;} else {index--;simp_h = - simp_h->parent;bdry_h=bdry_h->parent;} - } - */ - auto it_1 = simplices.simplex_vertex_range(s_1).begin(); auto end_1 = simplices.simplex_vertex_range(s_1).end(); - + // auto it_2 = simplices.simplex_vertex_range(s_2).begin(); auto end_2 = simplices.simplex_vertex_range(s_2).end(); - + // int orient = 0; int unmatch = 0; - + // for (int p = 0; p <= simplices.dimension(s_2); p++) { if (*it_1 != *it_2) { orient = p++; unmatch++; - } else { if (it_1 != end_1) { it_1++; @@ -269,6 +221,9 @@ class SimplicialComplex { return pow(-1, orient); } + /** + * @brief calculate the boundary matrices + */ void calculate_matrices() { boundary_matrices = std::vector(); for (int k = 0; k < simplices.dimension(); k++) { @@ -288,6 +243,4 @@ class SimplicialComplex { bool is_point(point_t pt) { return not(pt.size() == 0); } }; // class SimplicialComplex - -}; // namespace simplicial -#endif // SIMPLICIAL_COMPLEX +}; // namespace simplicial diff --git a/lib/scomplex/base_utils.hpp b/lib/scomplex/base_utils.hpp index 443666d..b538107 100644 --- a/lib/scomplex/base_utils.hpp +++ b/lib/scomplex/base_utils.hpp @@ -1,5 +1,4 @@ -#ifndef BASE_UTILS -#define BASE_UTILS +#pragma once #include #include @@ -18,7 +17,7 @@ std::list dedupe_list(std::list& list) { return no_reps_list; } -std::list dedupe_list(std::vector list) { +std::list dedupe_list(std::vector& list) { // O(n log(n)) std::set no_reps; std::list no_reps_list; @@ -31,7 +30,7 @@ std::list dedupe_list(std::vector list) { return no_reps_list; } -std::list dedupe_list(std::initializer_list list) { +std::list dedupe_list(std::initializer_list& list) { // O(n log(n)) std::set no_reps; std::list no_reps_list; @@ -43,5 +42,3 @@ std::list dedupe_list(std::initializer_list list) { } return no_reps_list; } - -#endif diff --git a/lib/scomplex/graph_utils.hpp b/lib/scomplex/graph_utils.hpp index 2672742..c5d408a 100644 --- a/lib/scomplex/graph_utils.hpp +++ b/lib/scomplex/graph_utils.hpp @@ -26,6 +26,7 @@ typedef typename boost::adjacency_list< // typedef typename boost::graph_traits::vertex_descriptor vertex_t; typedef std::vector point_t; +typedef std::vector simplex_t; Graph calculate_one_skelleton_graph( simplicial::SimplicialComplex& scomp // @@ -36,7 +37,7 @@ Graph calculate_one_skelleton_graph( std::vector> edges(num_edges); std::vector weights(num_edges); - for (std::vector edge : scomp.get_level(1)) { + for (simplex_t edge : scomp.get_level(1)) { int src_index = edge.at(0); int trg_index = edge.at(1); diff --git a/lib/scomplex/nn_utils.hpp b/lib/scomplex/nn_utils.hpp index b702403..756b4b3 100644 --- a/lib/scomplex/nn_utils.hpp +++ b/lib/scomplex/nn_utils.hpp @@ -4,18 +4,18 @@ #include #include - #include #include #include #include #include + #include typedef std::vector point_t; -// kd-tree point type -typedef CGAL::Cartesian_d::Point_d point_d; // d-dimensional point +// kd-tree point type, d-dimensional point +typedef CGAL::Cartesian_d::Point_d point_d; typedef CGAL::Search_traits_d> cartesian_traits; @@ -72,11 +72,13 @@ size_t nearest_neighbor_index(tree_t& tree, point_t point) { return ind; } -std::vector snap_points_to_indexes(tree_t& tree, - std::vector points) { +std::vector snap_points_to_indexes( // + tree_t& tree, // + std::vector points) { std::vector snapped_points; for (auto point : points) { - snapped_points.push_back(nearest_neighbor_index(tree, point)); + snapped_points.push_back( // + nearest_neighbor_index(tree, point)); } return snapped_points; } diff --git a/lib/scomplex/path_snapper.hpp b/lib/scomplex/path_snapper.hpp index 740aa4a..51693fe 100644 --- a/lib/scomplex/path_snapper.hpp +++ b/lib/scomplex/path_snapper.hpp @@ -11,6 +11,8 @@ class path_snapper { tree_t point_tree; // defined in nn_utils.hpp Graph vertex_graph; // defined in graph_utils.hpp + typedef std::vector, int>> geo_chain_t; + public: path_snapper(simplicial::SimplicialComplex& sc) { make_tree(point_tree, sc.points); @@ -27,6 +29,26 @@ class path_snapper { std::vector snapped_path = complete_path(vertex_graph, way_points); return snapped_path; - }; + } + + geo_chain_t get_chain_rep(std::vector path) { + auto vertex_path = snap_path(path); + geo_chain_t chain_rep; + for ( // + auto it = vertex_path.begin(); // + std::next(it) != vertex_path.end(); // + ++it) { + size_t s(*it), t(*(std::next(it))); + if (s < t) { + chain_rep.push_back( // + std::make_pair, int>( + std::vector{s, t}, 1)); + } else { + chain_rep.push_back( // + std::make_pair, int>( + std::vector{s, t}, 1)); + } + } + } }; }; From ad2e1035936d84febe9be1e5a45313c20cecb805 Mon Sep 17 00:00:00 2001 From: crvs Date: Tue, 7 Feb 2017 17:38:24 +0100 Subject: [PATCH 03/24] added vector support --- lib/scomplex/Simplicial_Complex.h | 31 ++++++------------------------- lib/scomplex/path_snapper.hpp | 25 ++++++++++++++++++------- src/libtest.cpp | 7 ------- 3 files changed, 24 insertions(+), 39 deletions(-) diff --git a/lib/scomplex/Simplicial_Complex.h b/lib/scomplex/Simplicial_Complex.h index ccfab6e..e768fa8 100644 --- a/lib/scomplex/Simplicial_Complex.h +++ b/lib/scomplex/Simplicial_Complex.h @@ -1,27 +1,3 @@ -/** - * @file Simplicial_Complex.h - * @brief utilities to create and manipulate geometric simplicial complexes - * @author Joao Carvalho - * @version 0.1 - * @date 2017-01-20 - * - * Copyright (C) - * 2017 - Joao Carvalho - * This program is free software; you can redistribute it and/or - * modify it under the terms of the GNU General Public License - * as published by the Free Software Foundation; either version 2 - * of the License, or (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - * - */ #pragma once #include @@ -182,6 +158,11 @@ class SimplicialComplex { */ bool is_geometric() { return geometric_q; } + size_t get_simplex_index(simplex_t simp) { + auto sh = simplices.find(simp); + return simplices.key(sh); + } + private: bool geometric_q; std::vector dimensions; @@ -196,7 +177,7 @@ class SimplicialComplex { return v_set; } - // need to do the boundary inclusion crap (you know what I mean) + // calculate the index of s_1 in the boundary of s_2 int boundary_index(Simplex_handle s_1, Simplex_handle s_2) { auto it_1 = simplices.simplex_vertex_range(s_1).begin(); auto end_1 = simplices.simplex_vertex_range(s_1).end(); diff --git a/lib/scomplex/path_snapper.hpp b/lib/scomplex/path_snapper.hpp index 51693fe..e4ea132 100644 --- a/lib/scomplex/path_snapper.hpp +++ b/lib/scomplex/path_snapper.hpp @@ -5,27 +5,27 @@ namespace snap { #include #include #include +#include + +typedef std::vector, int>> geo_chain_t; +typedef Eigen::SparseVector vector_t; class path_snapper { private: tree_t point_tree; // defined in nn_utils.hpp Graph vertex_graph; // defined in graph_utils.hpp - - typedef std::vector, int>> geo_chain_t; + simplicial::SimplicialComplex* s_comp; public: path_snapper(simplicial::SimplicialComplex& sc) { make_tree(point_tree, sc.points); vertex_graph = calculate_one_skelleton_graph(sc); + s_comp = ≻ } std::vector snap_path(std::vector path) { std::vector way_points = snap_points_to_indexes(point_tree, path); - std::copy(way_points.begin(), way_points.end(), - std::ostream_iterator(std::cout, " ")); - std::cout << '\n'; - std::vector snapped_path = complete_path(vertex_graph, way_points); return snapped_path; @@ -46,9 +46,20 @@ class path_snapper { } else { chain_rep.push_back( // std::make_pair, int>( - std::vector{s, t}, 1)); + std::vector{s, t}, -1)); } } } + + vector_t get_chain_vector(std::vector path) { + auto chain_rep = get_chain_rep(path); + vector_t vector_rep(s_comp->get_dimension(1)); + for (auto chain_el : chain_rep) { + simplex_t simp = std::get<0>(chain_el); + int val = std::get<1>(chain_el); + size_t index = s_comp->get_simplex_index(simp); + vector_rep.coeffRef(index) = val; + } + } }; }; diff --git a/src/libtest.cpp b/src/libtest.cpp index e3e5b64..673fec1 100644 --- a/src/libtest.cpp +++ b/src/libtest.cpp @@ -1,6 +1,5 @@ #include #include -//#include #include #include #include @@ -68,12 +67,6 @@ int main() { std::cout << "now quotient" << std::endl; auto sq = sc.quotient(f); - /* the points are working, what's up with the simplices?? - for (auto p : sq.points) { - std::cout << "point: " << p.at(0) << " " << p.at(1) << std::endl; - } - */ - range = sq.simplices.complex_simplex_range(); std::cout << "got this far" << std::endl; for (auto s : range) { From e28c7e3ed29e85d6522921ddca2c16a40306faab Mon Sep 17 00:00:00 2001 From: crvs Date: Wed, 8 Feb 2017 16:53:55 +0100 Subject: [PATCH 04/24] got rid of template argument --- lib/scomplex/Simplicial_Complex.h | 18 ++++++++-------- lib/scomplex/graph_utils.hpp | 5 ++--- lib/scomplex/path_snapper.hpp | 34 +++++++++++++++---------------- src/libtest.cpp | 4 ++-- 4 files changed, 30 insertions(+), 31 deletions(-) diff --git a/lib/scomplex/Simplicial_Complex.h b/lib/scomplex/Simplicial_Complex.h index e768fa8..1199b5f 100644 --- a/lib/scomplex/Simplicial_Complex.h +++ b/lib/scomplex/Simplicial_Complex.h @@ -11,9 +11,9 @@ #include #include +typedef std::vector point_t; namespace simplicial { -template class SimplicialComplex { private: // the type of matrix to be used @@ -52,7 +52,7 @@ class SimplicialComplex { int get_dimension(int level) { return dimensions.at(level); } - SimplicialComplex() { geometric_q = true; } + SimplicialComplex() { geometric_q = true; } /** * @brief builds a simplicial complex object out of points and "triangles" @@ -60,8 +60,8 @@ class SimplicialComplex { * @param points_a points * @param tris the top dimensional simplicial complex structure */ - SimplicialComplex(std::list points_a, - std::list> tris) { + SimplicialComplex(std::list points_a, + std::list> tris) { geometric_q = true; points = std::vector(); for (auto pt : points_a) { @@ -103,7 +103,7 @@ class SimplicialComplex { } } - SimplicialComplex quotient(int f(point_t)) { + SimplicialComplex quotient(int f(point_t)) { int n_points = 1; auto corresp = std::vector(); @@ -145,7 +145,7 @@ class SimplicialComplex { simp_list.push_back(dedupe_list(s_q)); } - auto quotient_sc = SimplicialComplex(q_points, simp_list); + auto quotient_sc = SimplicialComplex(q_points, simp_list); quotient_sc.geometric_q = geometric; return quotient_sc; @@ -169,10 +169,10 @@ class SimplicialComplex { std::vector boundary_matrices; - std::set vertex_set(Simplex_handle s) { - std::set v_set; + std::set vertex_set(Simplex_handle s) { + std::set v_set; for (auto v : simplices.simplex_vertex_range(s)) { - v_set.insert(s); + v_set.insert(v); } return v_set; } diff --git a/lib/scomplex/graph_utils.hpp b/lib/scomplex/graph_utils.hpp index c5d408a..7a2fc95 100644 --- a/lib/scomplex/graph_utils.hpp +++ b/lib/scomplex/graph_utils.hpp @@ -28,9 +28,8 @@ typedef typename boost::graph_traits::vertex_descriptor vertex_t; typedef std::vector point_t; typedef std::vector simplex_t; -Graph calculate_one_skelleton_graph( - simplicial::SimplicialComplex& scomp // - ) { +Graph calculate_one_skelleton_graph(simplicial::SimplicialComplex& scomp // + ) { int num_points{scomp.get_dimension(0)}; int num_edges{scomp.get_dimension(1)}; diff --git a/lib/scomplex/path_snapper.hpp b/lib/scomplex/path_snapper.hpp index e4ea132..fc68d6a 100644 --- a/lib/scomplex/path_snapper.hpp +++ b/lib/scomplex/path_snapper.hpp @@ -7,30 +7,15 @@ namespace snap { #include #include -typedef std::vector, int>> geo_chain_t; typedef Eigen::SparseVector vector_t; class path_snapper { private: tree_t point_tree; // defined in nn_utils.hpp Graph vertex_graph; // defined in graph_utils.hpp - simplicial::SimplicialComplex* s_comp; - - public: - path_snapper(simplicial::SimplicialComplex& sc) { - make_tree(point_tree, sc.points); - vertex_graph = calculate_one_skelleton_graph(sc); - s_comp = ≻ - } - - std::vector snap_path(std::vector path) { - std::vector way_points = - snap_points_to_indexes(point_tree, path); - std::vector snapped_path = - complete_path(vertex_graph, way_points); - return snapped_path; - } + simplicial::SimplicialComplex* s_comp; + typedef std::vector, int>> geo_chain_t; geo_chain_t get_chain_rep(std::vector path) { auto vertex_path = snap_path(path); geo_chain_t chain_rep; @@ -51,6 +36,21 @@ class path_snapper { } } + public: + path_snapper(simplicial::SimplicialComplex& sc) { + make_tree(point_tree, sc.points); + vertex_graph = calculate_one_skelleton_graph(sc); + s_comp = ≻ + } + + std::vector snap_path(std::vector path) { + std::vector way_points = + snap_points_to_indexes(point_tree, path); + std::vector snapped_path = + complete_path(vertex_graph, way_points); + return snapped_path; + } + vector_t get_chain_vector(std::vector path) { auto chain_rep = get_chain_rep(path); vector_t vector_rep(s_comp->get_dimension(1)); diff --git a/src/libtest.cpp b/src/libtest.cpp index 673fec1..ad88d91 100644 --- a/src/libtest.cpp +++ b/src/libtest.cpp @@ -30,7 +30,7 @@ int main() { simp_list.push_back(simp_t({0, 1, 2})); simp_list.push_back(simp_t({1, 2, 3})); - simplicial::SimplicialComplex sc(point_list, simp_list); + simplicial::SimplicialComplex sc(point_list, simp_list); for (auto pt : sc.points) { std::cout << pt.at(0) << " " << pt.at(1) << std::endl; } @@ -148,7 +148,7 @@ int main() { cells.push_back(cell); } - simplicial::SimplicialComplex rsc(points, cells); + simplicial::SimplicialComplex rsc(points, cells); snap::path_snapper snapper(rsc); From f689e9829fcd8b450d314ef49f7c152bcd74376b Mon Sep 17 00:00:00 2001 From: crvs Date: Thu, 9 Feb 2017 10:59:14 +0100 Subject: [PATCH 05/24] added direct constructor for path-snapper --- lib/scomplex/Simplicial_Complex.h | 10 +-- lib/scomplex/base_utils.hpp | 35 +++----- lib/scomplex/graph_utils.hpp | 3 - lib/scomplex/path_snapper.hpp | 17 +++- lib/scomplex/qhull_parsing.hpp | 13 +-- src/libtest.cpp | 134 ++---------------------------- 6 files changed, 47 insertions(+), 165 deletions(-) diff --git a/lib/scomplex/Simplicial_Complex.h b/lib/scomplex/Simplicial_Complex.h index 1199b5f..c0cb183 100644 --- a/lib/scomplex/Simplicial_Complex.h +++ b/lib/scomplex/Simplicial_Complex.h @@ -60,8 +60,8 @@ class SimplicialComplex { * @param points_a points * @param tris the top dimensional simplicial complex structure */ - SimplicialComplex(std::list points_a, - std::list> tris) { + SimplicialComplex(std::vector points_a, + std::vector> tris) { geometric_q = true; points = std::vector(); for (auto pt : points_a) { @@ -110,7 +110,7 @@ class SimplicialComplex { // making the correspondence between points in the original complex and // points in the quotient. point index 0 corresponds to the "virtual" // point (in case it exists). - std::list q_points; + std::vector q_points; bool geometric = true; for (auto p : points) { @@ -135,10 +135,10 @@ class SimplicialComplex { } // producing list of simplices - std::list> simp_list; + std::vector> simp_list; for (auto s : simplices.complex_simplex_range()) { - std::list s_q; + std::vector s_q; for (int v : simplices.simplex_vertex_range(s)) { s_q.push_back(corresp.at(v)); } diff --git a/lib/scomplex/base_utils.hpp b/lib/scomplex/base_utils.hpp index b538107..3867b7c 100644 --- a/lib/scomplex/base_utils.hpp +++ b/lib/scomplex/base_utils.hpp @@ -1,43 +1,30 @@ #pragma once +#include #include #include -#include -std::list dedupe_list(std::list& list) { +std::list dedupe_list(std::list& list) { // O(n log(n)) - std::set no_reps; - std::list no_reps_list; - for (int el : list) { + std::set no_reps; + std::list no_reps_list; + for (size_t el : list) { no_reps.insert(el); } - for (int el : no_reps) { + for (size_t el : no_reps) { no_reps_list.push_back(el); } return no_reps_list; } -std::list dedupe_list(std::vector& list) { +std::vector dedupe_list(std::vector& list) { // O(n log(n)) - std::set no_reps; - std::list no_reps_list; - for (int el : list) { + std::set no_reps; + std::vector no_reps_list; + for (size_t el : list) { no_reps.insert(el); } - for (int el : no_reps) { - no_reps_list.push_back(el); - } - return no_reps_list; -} - -std::list dedupe_list(std::initializer_list& list) { - // O(n log(n)) - std::set no_reps; - std::list no_reps_list; - for (int el : list) { - no_reps.insert(el); - } - for (int el : no_reps) { + for (size_t el : no_reps) { no_reps_list.push_back(el); } return no_reps_list; diff --git a/lib/scomplex/graph_utils.hpp b/lib/scomplex/graph_utils.hpp index 7a2fc95..e13346c 100644 --- a/lib/scomplex/graph_utils.hpp +++ b/lib/scomplex/graph_utils.hpp @@ -99,15 +99,12 @@ std::vector complete_path(const Graph& g, std::vector vec) { full_path.push_back(*s); while (t != vec.end()) { // get the next vertex - std::cout << "looking at edge" << *s << ", " << *t << "\n"; - // calculate and add the path between s and t (contains vertices other // than s) to the full path std::vector path_portion{shortest_path(g, *s, *t)}; for (vertex_t vert : path_portion) { full_path.push_back(vert); } - std::cout << "\n"; std::advance(s, 1); std::advance(t, 1); } diff --git a/lib/scomplex/path_snapper.hpp b/lib/scomplex/path_snapper.hpp index fc68d6a..266b0d2 100644 --- a/lib/scomplex/path_snapper.hpp +++ b/lib/scomplex/path_snapper.hpp @@ -14,6 +14,7 @@ class path_snapper { tree_t point_tree; // defined in nn_utils.hpp Graph vertex_graph; // defined in graph_utils.hpp simplicial::SimplicialComplex* s_comp; + bool owner; typedef std::vector, int>> geo_chain_t; geo_chain_t get_chain_rep(std::vector path) { @@ -38,9 +39,21 @@ class path_snapper { public: path_snapper(simplicial::SimplicialComplex& sc) { - make_tree(point_tree, sc.points); - vertex_graph = calculate_one_skelleton_graph(sc); + owner = false; s_comp = ≻ + make_tree(point_tree, s_comp->points); + vertex_graph = calculate_one_skelleton_graph(*s_comp); + } + + path_snapper(std::vector& pts, std::vector& cells) { + owner = true; + s_comp = new simplicial::SimplicialComplex(pts, cells); + make_tree(point_tree, s_comp->points); + vertex_graph = calculate_one_skelleton_graph(*s_comp); + } + + ~path_snapper() { + if (owner) delete s_comp; } std::vector snap_path(std::vector path) { diff --git a/lib/scomplex/qhull_parsing.hpp b/lib/scomplex/qhull_parsing.hpp index fdd384b..011c71b 100644 --- a/lib/scomplex/qhull_parsing.hpp +++ b/lib/scomplex/qhull_parsing.hpp @@ -11,7 +11,7 @@ #include typedef std::vector point_t; -typedef std::vector cell_t; +typedef std::vector cell_t; template std::vector tokenize(std::string str) { @@ -59,14 +59,14 @@ std::pair, std::vector> parse_qhull_file( std::getline(file_stream, line); // case of the first line of the file (get dimensions) if (line_number == 0) { - auto tokens = tokenize(line); + auto tokens = tokenize(line); dimensionality = tokens.at(0); // std::cout << "found first line! dimensions == " << dimensionality // << '\n'; } // second line of the file (get number of points) if (line_number == 1) { - auto tokens = tokenize(line); + auto tokens = tokenize(line); total_points = tokens.at(0); // std::cout << "found second line! number of points == " // << total_points << '\n'; @@ -79,7 +79,7 @@ std::pair, std::vector> parse_qhull_file( } // first line after points (get number of cells) if (line_number == total_points + 2) { - auto tokens = tokenize(line); + auto tokens = tokenize(line); total_cells = tokens.at(0); // std::cout << "found cell count! number of cells == " << // total_cells @@ -88,10 +88,11 @@ std::pair, std::vector> parse_qhull_file( // lines pertaining to cells if (line_number > total_points + 2 && line_number < total_points + total_cells + 3) { - auto cell = tokenize(line); + auto cell = tokenize(line); cells.push_back(cell); } } + file_stream.close(); return std::make_pair(points, cells); } @@ -109,7 +110,7 @@ int main() { } std::cout << cells.size() << '\n'; for (cell_t cell : cells) { - output_vector(cell); + output_vector(cell); } return 0; diff --git a/src/libtest.cpp b/src/libtest.cpp index ad88d91..bf5f61d 100644 --- a/src/libtest.cpp +++ b/src/libtest.cpp @@ -22,138 +22,22 @@ int main() { std::list simp_list; // --- - point_list.push_back(point_t({0, 0})); - point_list.push_back(point_t({0, 5})); - point_list.push_back(point_t({1, 0})); - point_list.push_back(point_t({1, 1})); - - simp_list.push_back(simp_t({0, 1, 2})); - simp_list.push_back(simp_t({1, 2, 3})); - - simplicial::SimplicialComplex sc(point_list, simp_list); - for (auto pt : sc.points) { - std::cout << pt.at(0) << " " << pt.at(1) << std::endl; - } - - auto range = sc.simplices.complex_simplex_range(); - for (auto s : range) { - std::cout << sc.simplices.dimension(s) << " , " << sc.simplices.key(s) - << " : "; - for (auto v : sc.simplices.simplex_vertex_range(s)) { - std::cout << v << " "; - } - // iterate over boundary - std::cout << std::endl - << " boundary: "; - for (auto sb : sc.simplices.boundary_simplex_range(s)) { - auto v_range = sc.simplices.simplex_vertex_range(sb); - for (auto v : v_range) { - std::cout << v << " "; - } - std::cout << " ; "; - } - std::cout << std::endl - << " coboundary: "; - for (auto sb : sc.simplices.cofaces_simplex_range(s, 1)) { - auto v_range = sc.simplices.simplex_vertex_range(sb); - for (auto v : v_range) { - std::cout << v << " "; - } - std::cout << " ; "; - } - std::cout << std::endl; - } - - std::cout << "now quotient" << std::endl; - auto sq = sc.quotient(f); - - range = sq.simplices.complex_simplex_range(); - std::cout << "got this far" << std::endl; - for (auto s : range) { - std::cout << sq.simplices.dimension(s) << " , " << sq.simplices.key(s) - << " : "; - for (auto v : sq.simplices.simplex_vertex_range(s)) { - std::cout << v << " "; - } - // iterate over boundary - std::cout << std::endl - << " boundary: "; - for (auto sb : sq.simplices.boundary_simplex_range(s)) { - auto v_range = sq.simplices.simplex_vertex_range(sb); - for (auto v : v_range) { - std::cout << v << " "; - } - std::cout << " ; "; - } - std::cout << std::endl - << " coboundary: "; - for (auto sb : sq.simplices.cofaces_simplex_range(s, 1)) { - auto v_range = sq.simplices.simplex_vertex_range(sb); - for (auto v : v_range) { - std::cout << v << " "; - } - std::cout << " ; "; - } - std::cout << std::endl; - } - std::cout << "got this far" << std::endl; - - // checking the matrices - auto g = sc.get_boundary(1); - std::cout << g << std::endl; - g = sc.get_boundary(0); - std::cout << g << std::endl; - - g = sq.get_boundary(1); - std::cout << g << std::endl; - g = sq.get_boundary(0); - std::cout << g << std::endl; - - std::cout << "printing level: " << std::endl; - for (auto s : sq.get_level(1)) { - std::cout << s.size() << ": "; - std::cout << s.at(0) << " "; - std::cout << s.at(1) << " "; - std::cout << std::endl; - } - - auto G = calculate_one_skelleton_graph(sc); - auto p = shortest_path(G, 0, 3); - - typedef typename decltype(G)::vertex_descriptor vertex_t; - - std::cout << '\n'; - std::copy(p.begin(), p.end(), - std::ostream_iterator{std::cout, " "}); - std::cout << '\n'; - - std::vector point_vec{point_list.begin(), point_list.end()}; - tree_t rt; - make_tree(rt, point_vec); - - point_t q{0.1, 0.2}; - point_t q_prime = nearest_neighbor(rt, q); - // get_r_tree_from_points<2>(point_list); - std::ostream_iterator outstr(std::cout, " "); - std::copy(q_prime.begin(), q_prime.end(), outstr); - std::vector points_v; std::vector cells_v; - std::tie(points_v, cells_v) = - parse_qhull_file("/home/crvs/dev-cpp/examples/qhull-test/qh-test.dat"); - std::list points(points_v.begin(), points_v.end()); - std::list> cells; - for (auto cell_v : cells_v) { - std::list cell(cell_v.begin(), cell_v.end()); - cells.push_back(cell); - } - simplicial::SimplicialComplex rsc(points, cells); + std::string filename; + std::cout << "insert file name: "; + std::cin >> filename; + std::tie(points_v, cells_v) = + // parse_qhull_file("/home/crvs/dev-cpp/examples/qhull-test/qh-test.dat"); + parse_qhull_file(filename); + simplicial::SimplicialComplex rsc(points_v, cells_v); snap::path_snapper snapper(rsc); + // snap::path_snapper snapper(points_v, cells_v); std::cout << "\n"; - auto path = snapper.snap_path({{-1, -1}, {1, 1}}); + auto path = snapper.snap_path({{-1, -1}, {1, 0}, {1, 1}}); std::copy(path.begin(), path.end(), std::ostream_iterator(std::cout, " ")); return 0; From b0fabd6ef017b9a4116ccd168fea796734e888f1 Mon Sep 17 00:00:00 2001 From: crvs Date: Tue, 14 Feb 2017 14:41:53 +0100 Subject: [PATCH 06/24] added chain calculations --- CMakeLists.txt | 2 +- lib/scomplex/Simplicial_Complex.h | 36 +++--- lib/scomplex/chain-calc.hpp | 181 ++++++++++++++++++++++++++++++ lib/scomplex/path_snapper.hpp | 2 + lib/scomplex/qhull_parsing.hpp | 46 +------- src/libtest.cpp | 40 +++++-- 6 files changed, 241 insertions(+), 66 deletions(-) create mode 100644 lib/scomplex/chain-calc.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 7eeadf1..f3f4fac 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,7 +10,7 @@ project(simplicial CXX) # * boost under /usr/include/boost # * CGAL under /usr/include/CGAL # -set(EIGEN3_INCLUDE_DIR "/usr/include/eigen3") +set(EIGEN3_INCLUDE_DIR "/usr/local/include/eigen3") include_directories(${EIGEN3_INCLUDE_DIR}) set(CGAL_INCLUDE_DIR "/usr/local/include") include_directories(${CGAL_INCLUDE_DIR}) diff --git a/lib/scomplex/Simplicial_Complex.h b/lib/scomplex/Simplicial_Complex.h index c0cb183..17a683b 100644 --- a/lib/scomplex/Simplicial_Complex.h +++ b/lib/scomplex/Simplicial_Complex.h @@ -14,11 +14,13 @@ typedef std::vector point_t; namespace simplicial { +class NoChain {}; + class SimplicialComplex { private: // the type of matrix to be used typedef typename Eigen::SparseMatrix matrix_t; - typedef typename Eigen::SparseVector chain_t; + typedef typename Eigen::SparseVector vector_t; // simplex tree options struct SimpleOptions : Gudhi::Simplex_tree_options_full_featured { @@ -52,7 +54,7 @@ class SimplicialComplex { int get_dimension(int level) { return dimensions.at(level); } - SimplicialComplex() { geometric_q = true; } + SimplicialComplex() : geometric_q(true) {} /** * @brief builds a simplicial complex object out of points and "triangles" @@ -61,12 +63,8 @@ class SimplicialComplex { * @param tris the top dimensional simplicial complex structure */ SimplicialComplex(std::vector points_a, - std::vector> tris) { - geometric_q = true; - points = std::vector(); - for (auto pt : points_a) { - points.push_back(pt); - } + std::vector> tris) + : geometric_q(true), points(points_a) { for (auto s : tris) { // cant have simplices with repeated vertices so we dedupe the lists simplices.insert_simplex_and_subfaces(dedupe_list(s)); @@ -93,13 +91,16 @@ class SimplicialComplex { } matrix_t get_boundary(int d) { + // uninstantiated boundary matrices if (boundary_matrices.size() == 0) { calculate_matrices(); } - if (d < boundary_matrices.size()) { + + // now they have to be instantiated, get them + if (0 <= d && d < boundary_matrices.size()) { return boundary_matrices.at(d); } else { - return matrix_t(0, 0); + throw NoChain(); } } @@ -107,8 +108,10 @@ class SimplicialComplex { int n_points = 1; auto corresp = std::vector(); - // making the correspondence between points in the original complex and - // points in the quotient. point index 0 corresponds to the "virtual" + // making the correspondence between points in the original complex + // and + // points in the quotient. point index 0 corresponds to the + // "virtual" // point (in case it exists). std::vector q_points; bool geometric = true; @@ -119,7 +122,8 @@ class SimplicialComplex { q_points.push_back(p); } else { if (geometric) { - // need to realize that it is not geometric, (flip geometric + // need to realize that it is not geometric, (flip + // geometric // and add the virtual point auto pos = q_points.begin(); q_points.insert(pos, point_t()); @@ -163,6 +167,8 @@ class SimplicialComplex { return simplices.key(sh); } + int dimension() { return simplices.dimension(); } + private: bool geometric_q; std::vector dimensions; @@ -221,7 +227,9 @@ class SimplicialComplex { } } - bool is_point(point_t pt) { return not(pt.size() == 0); } + bool is_point(point_t pt) { + return (geometric_q ? true : not(pt.size() == 0)); + } }; // class SimplicialComplex }; // namespace simplicial diff --git a/lib/scomplex/chain-calc.hpp b/lib/scomplex/chain-calc.hpp new file mode 100644 index 0000000..0e4dc47 --- /dev/null +++ b/lib/scomplex/chain-calc.hpp @@ -0,0 +1,181 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +typedef Eigen::SparseVector vector_t; +typedef std::pair chain_t; +typedef Eigen::SparseMatrix matrix_t; +// typedef std::vector cell_t; + +class bounding_chain { + std::vector boundary_matrices; + simplicial::SimplicialComplex* s_comp; + bool owner; + void populate_matrices(); + + public: + bounding_chain(simplicial::SimplicialComplex& sc); + bounding_chain(std::vector& points, std::vector& tris); + ~bounding_chain(); + chain_t get_bounding_chain(chain_t&); + void draw_chain(std::string, chain_t); +}; + +class NonZeroChain : public std::exception {}; + +//--------------------------------------- +// implementation + +bounding_chain::~bounding_chain() { + if (this->owner) delete this->s_comp; + + // we always own the matrices + for (int i = this->boundary_matrices.size() - 1; i >= 0; --i) + delete this->boundary_matrices.at(i); +} + +bounding_chain::bounding_chain(simplicial::SimplicialComplex& sc) { + this->owner = false; + this->s_comp = ≻ + this->populate_matrices(); +} + +bounding_chain::bounding_chain(std::vector& points, + std::vector& tris) { + this->owner = true; + this->s_comp = new simplicial::SimplicialComplex(points, tris); + this->populate_matrices(); +} + +void bounding_chain::populate_matrices() { + for (int d = 0; d < this->s_comp->dimension(); ++d) { + auto my_matrix = this->s_comp->get_boundary(d); + this->boundary_matrices.push_back(new matrix_t(my_matrix)); + matrix_t* put_matrix = this->boundary_matrices.at(d); + } +} +// +// round vectors for comparison +vector_t round_vec(vector_t vec) { + vector_t rounded_vec(vec.rows()); + for (int i = 0; i < vec.rows(); ++i) { + int coef = round(vec.coeffRef(i)); + if (coef != 0) rounded_vec.coeffRef(i) = coef; + } + return rounded_vec; +} + +bool equals(vector_t vec1, vector_t vec2) { + for (int i = 0; i < vec1.rows(); ++i) { + if (vec1.coeffRef(i) != vec2.coeffRef(i)) { + std::cout << i << ": " << vec1.coeffRef(i) << " " + << vec2.coeffRef(i) << '\n'; + return false; + } + } + return true; +} + +void bounding_chain::draw_chain(std::string filename, chain_t chain) { + typedef boost::geometry::model::d2::point_xy point_type; + typedef boost::geometry::model::ring polygon_type; + + std::vector polys; + + for (auto tri : this->s_comp->get_level(2)) { + std::vector polygon_vec; + for (auto k : tri) { + point_t pt = this->s_comp->points.at(k); + point_type pt_g(pt[0], pt[1]); + polygon_vec.push_back(pt_g); + } + polygon_type poly(polygon_vec.begin(), polygon_vec.end()); + polys.push_back(poly); + } + + // now we want to draw this + std::ofstream svg(filename.c_str()); + + boost::geometry::svg_mapper mapper(svg, 150, 150); + + int chain_d = std::get<0>(chain); + vector_t chain_v = std::get<1>(chain); + for (int i = 0; i < chain_v.rows(); ++i) { + polygon_type poly = polys.at(i); + mapper.add(poly); + if (chain_d == 2) { + if (chain_v.coeffRef(i) == 0) + mapper.map( + poly, + "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153," + "204,0);stroke-width:2"); + else + mapper.map( + poly, + "fill-opacity:0.5;fill:rgb(230,104,58);stroke:rgb(230," + "104,58);stroke-width:2"); + } else { + mapper.map(poly, + "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153," + "204,0);stroke-width:2"); + } + + if (chain_d == 1) { + auto lines = this->s_comp->get_level(1); + for (int i = 0; i < lines.size(); ++i) { + if (chain_v.coeffRef(i) != 0) { + std::vector line_vec; + for (auto k : lines.at(i)) { + point_t pt = this->s_comp->points.at(k); + point_type pt_g(pt[0], pt[1]); + line_vec.push_back(pt_g); + } + polygon_type line_g(line_vec.begin(), line_vec.end()); + mapper.map( + line_g, + "fill-opacity:0.5;fill:rgb(250,0,0);stroke:rgb(250," + "0,0);stroke-width:2"); + } + } + } + } +} + +chain_t bounding_chain::get_bounding_chain(chain_t& chain) { + int chain_d; + vector_t chain_v; + std::tie(chain_d, chain_v) = chain; + + if (chain_d >= this->s_comp->dimension()) throw NonZeroChain(); + + Eigen::LeastSquaresConjugateGradient lscg; + lscg.compute(*(this->boundary_matrices.at(chain_d))); + + vector_t bound_chain(round_vec(lscg.solve(chain_v))); + + vector_t result( + round_vec(*(this->boundary_matrices.at(chain_d)) * bound_chain)); + if (equals(result, chain_v)) + return chain_t(chain_d - 1, bound_chain); + else { + std::cout << "where to write the svg file?"; + std::string filename; + std::cin >> filename; + this->draw_chain(filename, chain_t(1, bound_chain)); + std::cout << '\n'; + std::cout << result << '\n'; + std::cout << chain_v << '\n'; + // throw NonZeroChain(); + } +} diff --git a/lib/scomplex/path_snapper.hpp b/lib/scomplex/path_snapper.hpp index 266b0d2..ca1d9b2 100644 --- a/lib/scomplex/path_snapper.hpp +++ b/lib/scomplex/path_snapper.hpp @@ -35,6 +35,7 @@ class path_snapper { std::vector{s, t}, -1)); } } + return chain_rep; } public: @@ -73,6 +74,7 @@ class path_snapper { size_t index = s_comp->get_simplex_index(simp); vector_rep.coeffRef(index) = val; } + return vector_rep; } }; }; diff --git a/lib/scomplex/qhull_parsing.hpp b/lib/scomplex/qhull_parsing.hpp index 011c71b..d9e5b13 100644 --- a/lib/scomplex/qhull_parsing.hpp +++ b/lib/scomplex/qhull_parsing.hpp @@ -16,7 +16,8 @@ typedef std::vector cell_t; template std::vector tokenize(std::string str) { std::istringstream sstream(str); - std::vector tokens{std::istream_iterator{sstream}, + std::vector tokens{// + std::istream_iterator{sstream}, std::istream_iterator{}}; return tokens; } @@ -38,18 +39,11 @@ std::pair, std::vector> parse_qhull_file( std::string line; - int current_line_number{0}; - int dimensionality; - // - int num_points; - int num_cells; - // - int cur_point_num; - int cur_cell_num; + int current_line_number{0}, dimensionality; + int num_points, num_cells; + int cur_point_num, cur_cell_num; - int line_number{0}; - int total_points{0}; - int total_cells{0}; + int line_number{0}, total_points{0}, total_cells{0}; std::vector points; std::vector cells; @@ -61,15 +55,11 @@ std::pair, std::vector> parse_qhull_file( if (line_number == 0) { auto tokens = tokenize(line); dimensionality = tokens.at(0); - // std::cout << "found first line! dimensions == " << dimensionality - // << '\n'; } // second line of the file (get number of points) if (line_number == 1) { auto tokens = tokenize(line); total_points = tokens.at(0); - // std::cout << "found second line! number of points == " - // << total_points << '\n'; } // lines pertaining to points if (line_number > 1 && line_number < total_points + 2) { @@ -81,9 +71,6 @@ std::pair, std::vector> parse_qhull_file( if (line_number == total_points + 2) { auto tokens = tokenize(line); total_cells = tokens.at(0); - // std::cout << "found cell count! number of cells == " << - // total_cells - // << '\n'; } // lines pertaining to cells if (line_number > total_points + 2 && @@ -95,24 +82,3 @@ std::pair, std::vector> parse_qhull_file( file_stream.close(); return std::make_pair(points, cells); } - -/* - * usage example - * -int main() { - std::vector points; - std::vector cells; - std::tie(points, cells) = parse_qhull_file("qh-test.dat"); - - std::cout << points.size() << '\n'; - for (point_t point : points) { - output_vector(point); - } - std::cout << cells.size() << '\n'; - for (cell_t cell : cells) { - output_vector(cell); - } - - return 0; -} -*/ diff --git a/src/libtest.cpp b/src/libtest.cpp index bf5f61d..c9ca754 100644 --- a/src/libtest.cpp +++ b/src/libtest.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include "boost/graph/graph_traits.hpp" @@ -12,33 +13,50 @@ #include typedef std::vector point_t; -typedef std::list simp_t; -typedef std::pair> chain_t; +typedef std::vector cell_t; int f(point_t a) { return a.at(0) < .5 ? 0 : 1; } int main() { - std::list point_list; - std::list simp_list; - // --- - std::vector points_v; std::vector cells_v; std::string filename; std::cout << "insert file name: "; std::cin >> filename; - std::tie(points_v, cells_v) = - // parse_qhull_file("/home/crvs/dev-cpp/examples/qhull-test/qh-test.dat"); - parse_qhull_file(filename); + std::tie(points_v, cells_v) = parse_qhull_file(filename); simplicial::SimplicialComplex rsc(points_v, cells_v); + auto qsc(rsc.quotient(f)); snap::path_snapper snapper(rsc); - // snap::path_snapper snapper(points_v, cells_v); - std::cout << "\n"; auto path = snapper.snap_path({{-1, -1}, {1, 0}, {1, 1}}); std::copy(path.begin(), path.end(), std::ostream_iterator(std::cout, " ")); + std::cout << '\n'; + + bounding_chain chain_calc(points_v, cells_v); + + path = snapper.snap_path({{-1, -1}, {1, 0}, {1, 1}, {-1, -1}}); + std::copy(path.begin(), path.end(), + std::ostream_iterator(std::cout, " ")); + std::cout << '\n'; + vector_t path_v( + snapper.get_chain_vector({{-1, -1}, {1, 0}, {1, 1}, {-1, -1}})); + + // try { + chain_t chain(2, path_v); + chain = chain_calc.get_bounding_chain(chain); + std::cout << path_v << '\n'; + vector_t b_chain = round_vec(std::get<1>(chain)); + std::cout << b_chain; + // std::cout << chain_r << '\n'; + std::cout << round_vec(rsc.get_boundary(1) * std::get<1>(chain)) << '\n'; + std::cout << round_vec(rsc.get_boundary(1) * b_chain) << '\n'; + + //} catch (NonZeroChain) { + std::cout << "not a zero chain"; + //} + return 0; } From 29db7b1e0668e876811cae7da1a29ac98a9eafde Mon Sep 17 00:00:00 2001 From: crvs Date: Tue, 14 Feb 2017 18:46:51 +0100 Subject: [PATCH 07/24] refactored simplicial complex --- CMakeLists.txt | 20 +-- lib/scomplex/base_utils.hpp | 10 +- lib/scomplex/simplicial_complex.cpp | 222 ++++++++++++++++++++++++++++ lib/scomplex/simplicial_complex.hpp | 35 +++++ lib/scomplex/types.hpp | 15 ++ src/impl_test.cpp | 19 +++ 6 files changed, 310 insertions(+), 11 deletions(-) create mode 100644 lib/scomplex/simplicial_complex.cpp create mode 100644 lib/scomplex/simplicial_complex.hpp create mode 100644 lib/scomplex/types.hpp create mode 100644 src/impl_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index f3f4fac..e54d9c3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -23,19 +23,21 @@ include(${CGAL_USE_FILE}) set(CMAKE_BUILD_TYPE Debug) -include_directories("./lib") +include_directories("./lib") # where the header files are set(CMAKE_CXX_FLAGS "-std=c++11 ${CMAKE_CXX_FLAGS}") +set(CMAKE_CXX_FLAGS "-lspatialindex_c -lspatialindex ${CMAKE_CXX_FLAGS}") +add_library(test SHARED "./lib/scomplex/simplicial_complex.cpp") -set(STUPID_SOURCES "./src/main.cpp") -add_executable(my_stupid_test ${STUPID_SOURCES}) +add_executable(impl_test "src/impl_test.cpp") +target_link_libraries(impl_test test) -set(YOTHER_STUPID_SOURCES "./src/graphtest.cpp") -add_executable(yother_stupid_test ${YOTHER_STUPID_SOURCES}) +# set(STUPID_SOURCES "./src/main.cpp") +# add_executable(my_stupid_test ${STUPID_SOURCES}) -# set(CMAKE_CXX_FLAGS "-lspatialindex_c -lspatialindex ${CMAKE_CXX_FLAGS}") -set(OTHER_STUPID_SOURCES "./src/libtest.cpp") -add_executable(my_other_stupid_test ${OTHER_STUPID_SOURCES}) +# set(YOTHER_STUPID_SOURCES "./src/graphtest.cpp") +# add_executable(yother_stupid_test ${YOTHER_STUPID_SOURCES}) -add_executable(nn_test "./src/nn_test") +# set(OTHER_STUPID_SOURCES "./src/libtest.cpp") +# add_executable(my_other_stupid_test ${OTHER_STUPID_SOURCES}) diff --git a/lib/scomplex/base_utils.hpp b/lib/scomplex/base_utils.hpp index 3867b7c..e5471fe 100644 --- a/lib/scomplex/base_utils.hpp +++ b/lib/scomplex/base_utils.hpp @@ -17,11 +17,11 @@ std::list dedupe_list(std::list& list) { return no_reps_list; } -std::vector dedupe_list(std::vector& list) { +std::vector dedupe_vec(std::vector& vec) { // O(n log(n)) std::set no_reps; std::vector no_reps_list; - for (size_t el : list) { + for (size_t el : vec) { no_reps.insert(el); } for (size_t el : no_reps) { @@ -29,3 +29,9 @@ std::vector dedupe_list(std::vector& list) { } return no_reps_list; } + +template +std::unique_ptr new_unique(T* t) { + T* tt = new T(t); + return std::unique_ptr(tt); +} diff --git a/lib/scomplex/simplicial_complex.cpp b/lib/scomplex/simplicial_complex.cpp new file mode 100644 index 0000000..087e38f --- /dev/null +++ b/lib/scomplex/simplicial_complex.cpp @@ -0,0 +1,222 @@ +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace gsimp { + +struct simplicial_complex::impl { + // auxiliary types + struct SimpleOptions : Gudhi::Simplex_tree_options_full_featured { + typedef size_t Vertex_handle; + }; + typedef Gudhi::Simplex_tree simp_tree; + typedef Gudhi::Simplex_tree::Simplex_handle simp_handle; + typedef std::vector level_t; + typedef std::vector> levels_t; + + // member variables + std::vector points; + simp_tree simplices; + bool quotient_q; + std::vector boundary_matrices; + levels_t levels; + + impl() : quotient_q(false){}; + + impl(std::vector& arg_points, std::vector& arg_tris) + : quotient_q(false), points(arg_points) { + // create the simplex tree + for (auto tri : arg_tris) { + simplices.insert_simplex_and_subfaces(dedupe_vec(tri)); + int d = tri.size() - 1; + if (simplices.dimension() < d) simplices.set_dimension(d); + } + + // assign a key to each simplex in each level + std::vector count(simplices.dimension() + 1, 0); + for (int i = 0; i < simplices.dimension() + 1; ++i) { + level_t* level = new level_t(); + levels.push_back(std::unique_ptr(level)); + } + + auto simplex_range = simplices.complex_simplex_range(); + for (auto s : simplex_range) { + int d = simplices.dimension(s); + simplices.assign_key(s, count.at(d)++); + levels.at(d)->push_back(&s); + } + } + + ~impl() { + // get rid of the levels so they won't dangle + for (int i = levels.size(); i <= 0; --i) levels.at(i).reset(); + }; + + size_t get_level_size(int level) { return levels.at(level)->size(); } + + // calculate the index of s_1 in the boundary of s_2 + int boundary_index(simp_handle s_1, simp_handle s_2) { + auto it_1 = simplices.simplex_vertex_range(s_1).begin(); + auto end_1 = simplices.simplex_vertex_range(s_1).end(); + // + auto it_2 = simplices.simplex_vertex_range(s_2).begin(); + auto end_2 = simplices.simplex_vertex_range(s_2).end(); + // + int orient = 0; + int unmatch = 0; + // + for (int p = 0; p <= simplices.dimension(s_2); p++) { + if (*it_1 != *it_2) { + orient = p++; + unmatch++; + } else { + if (it_1 != end_1) it_1++; + it_2++; + } + } + return pow(-1, orient); + } + + void calculate_matrices() { + boundary_matrices = std::vector(); + for (int k = 0; k < simplices.dimension(); k++) { + boundary_matrices.push_back( + matrix_t(get_level_size(k), get_level_size(k + 1))); + } + for (auto s : simplices.complex_simplex_range()) { + int j = simplices.key(s); + for (auto bs : simplices.boundary_simplex_range(s)) { + int i = simplices.key(bs); + int k = simplices.dimension(bs); + boundary_matrices.at(k).coeffRef(i, j) = boundary_index(bs, s); + } + } + } + + std::vector get_level(int level) { + std::vector level_cells; + for (auto simp : *levels.at(level)) { + cell_t v_simp; + for (auto v : simplices.simplex_vertex_range(*simp)) { + v_simp.push_back(v); + } + level_cells.push_back(v_simp); + } + return level_cells; + } +}; // struct impl + +int simplicial_complex::get_level_size(int level) { + return p_impl->get_level_size(level); +} + +simplicial_complex::simplicial_complex(std::vector& arg_points, + std::vector& arg_tris) { + p_impl = std::shared_ptr(new impl(arg_points, arg_tris)); +} + +simplicial_complex::~simplicial_complex() {} + +std::vector simplicial_complex::get_points() { return p_impl->points; } + +std::vector simplicial_complex::get_level(int level) { + return p_impl->get_level(level); +} + +matrix_t simplicial_complex::get_boundary_matrix(int d) { + // uninstantiated boundary matrices + if (p_impl->boundary_matrices.size() == 0) { + p_impl->calculate_matrices(); + } + + // now they have to be instantiated, get them + if (0 <= d && d < p_impl->boundary_matrices.size()) { + return p_impl->boundary_matrices.at(d); + } else { + throw NoChain(); + } +} + +int simplicial_complex::dimension() { return p_impl->simplices.dimension(); } + +size_t simplicial_complex::get_index_of_simplex(cell_t simp) { + auto sh = p_impl->simplices.find(simp); + return p_impl->simplices.key(sh); +} + +bool simplicial_complex::is_quotient() { return p_impl->quotient_q; } + +simplicial_complex simplicial_complex::quotient(int char_fun(point_t)) { + int n_points = 1; + auto corresp = std::vector(); + + // making the correspondence between points in the original complex + // and points in the quotient. point index 0 corresponds to the + // "virtual" point (in case it exists). + std::vector points; + bool quotient_q = p_impl->quotient_q; + + for (auto p : p_impl->points) { + if (char_fun(p) != 0) { + corresp.push_back(n_points++); + points.push_back(p); + } else { + if (not quotient_q) { + // need to realize that it is not geometric, (flip + // quotient_q and add the virtual point + points.insert(points.begin(), point_t()); + quotient_q = false; + } + corresp.push_back(0); + } + } + + if (not quotient_q) { + for (int i = 0; i < corresp.size(); i++) { + corresp.at(i) -= 1; + } + } + + // producing list of simplices + std::vector simp_list; + + for (auto s : p_impl->simplices.complex_simplex_range()) { + cell_t s_q; + for (int v : p_impl->simplices.simplex_vertex_range(s)) + s_q.push_back(corresp.at(v)); + simp_list.push_back(s_q); + } + + auto quotient_sc = simplicial_complex(points, simp_list); + quotient_sc.p_impl->quotient_q = quotient_q; + + return quotient_sc; +} + +}; // namespace gsimp + +// DEPRECATED + +/* + *std::set vertex_set(simp_handle s) { + * std::set v_set; + * for (auto v : simplices.simplex_vertex_range(s)) { + * v_set.insert(v); + * } + * return v_set; + *} + */ + +/* MAY BE INCLUDED AGAIN WHEN WE START DOING MORE WITH QUOTIENTS + * + * bool is_point(point_t pt) { + * return (not quotient_q ? true : not(pt.size() == 0)); + * } + * + */ diff --git a/lib/scomplex/simplicial_complex.hpp b/lib/scomplex/simplicial_complex.hpp new file mode 100644 index 0000000..14c20a9 --- /dev/null +++ b/lib/scomplex/simplicial_complex.hpp @@ -0,0 +1,35 @@ +#pragma once + +#include +#include + +namespace gsimp { + +using namespace gsimp; + +class NoChain {}; + +class simplicial_complex { + private: + // implementation details + struct impl; + std::shared_ptr p_impl; + + public: + // constructor (no default) + simplicial_complex(std::vector&, std::vector&); + + // destructor + ~simplicial_complex(); + + std::vector get_points(); + std::vector get_level(int); + int dimension(); + int get_level_size(int); + matrix_t get_boundary_matrix(int); + bool is_quotient(); + size_t get_index_of_simplex(cell_t); + + simplicial_complex quotient(int char_fun(point_t)); +}; // class simplicial_complex +}; // namespace gsimp diff --git a/lib/scomplex/types.hpp b/lib/scomplex/types.hpp new file mode 100644 index 0000000..3cf2783 --- /dev/null +++ b/lib/scomplex/types.hpp @@ -0,0 +1,15 @@ +#include +#include + +namespace gsimp { +// geometric types +typedef typename std::vector point_t; +typedef typename std::vector cell_t; + +// linear algebra types +typedef typename Eigen::SparseMatrix matrix_t; +typedef typename Eigen::SparseVector vector_t; + +// chains +typedef typename std::pair chain_t; +}; diff --git a/src/impl_test.cpp b/src/impl_test.cpp new file mode 100644 index 0000000..a7403d2 --- /dev/null +++ b/src/impl_test.cpp @@ -0,0 +1,19 @@ +#include + +#include +#include + +int main() { + std::vector points_v; + std::vector cells_v; + + std::string filename; + std::cout << "insert file name: "; + std::cin >> filename; + std::tie(points_v, cells_v) = parse_qhull_file(filename); + + gsimp::simplicial_complex s_comp(points_v, cells_v); + std::cout << s_comp.dimension(); + + return 0; +} From d37cd2547e67bd709035c537866d2e382054a902 Mon Sep 17 00:00:00 2001 From: crvs Date: Tue, 14 Feb 2017 20:10:08 +0100 Subject: [PATCH 08/24] fixed the quotient --- lib/scomplex/simplicial_complex.cpp | 33 +++++++++++++++++------------ lib/scomplex/simplicial_complex.hpp | 8 +++---- lib/scomplex/types.hpp | 4 ++++ src/impl_test.cpp | 17 +++++++++++---- 4 files changed, 41 insertions(+), 21 deletions(-) diff --git a/lib/scomplex/simplicial_complex.cpp b/lib/scomplex/simplicial_complex.cpp index 087e38f..98813f5 100644 --- a/lib/scomplex/simplicial_complex.cpp +++ b/lib/scomplex/simplicial_complex.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include namespace gsimp { @@ -121,6 +122,16 @@ simplicial_complex::simplicial_complex(std::vector& arg_points, p_impl = std::shared_ptr(new impl(arg_points, arg_tris)); } +simplicial_complex::simplicial_complex(const simplicial_complex& other) { + p_impl = other.p_impl; +} + +simplicial_complex& simplicial_complex::operator=( + const simplicial_complex& other) { + p_impl = other.p_impl; + return *this; +} + simplicial_complex::~simplicial_complex() {} std::vector simplicial_complex::get_points() { return p_impl->points; } @@ -131,16 +142,13 @@ std::vector simplicial_complex::get_level(int level) { matrix_t simplicial_complex::get_boundary_matrix(int d) { // uninstantiated boundary matrices - if (p_impl->boundary_matrices.size() == 0) { - p_impl->calculate_matrices(); - } + if (p_impl->boundary_matrices.size() == 0) p_impl->calculate_matrices(); // now they have to be instantiated, get them - if (0 <= d && d < p_impl->boundary_matrices.size()) { + if (0 <= d && d < p_impl->boundary_matrices.size()) return p_impl->boundary_matrices.at(d); - } else { + else throw NoChain(); - } } int simplicial_complex::dimension() { return p_impl->simplices.dimension(); } @@ -152,6 +160,8 @@ size_t simplicial_complex::get_index_of_simplex(cell_t simp) { bool simplicial_complex::is_quotient() { return p_impl->quotient_q; } +// TODO: get top dimensional cells from the simplicial complex to add to the +// simplex tree (instead of all the simplices). simplicial_complex simplicial_complex::quotient(int char_fun(point_t)) { int n_points = 1; auto corresp = std::vector(); @@ -171,17 +181,14 @@ simplicial_complex simplicial_complex::quotient(int char_fun(point_t)) { // need to realize that it is not geometric, (flip // quotient_q and add the virtual point points.insert(points.begin(), point_t()); - quotient_q = false; + quotient_q = true; } corresp.push_back(0); } } - if (not quotient_q) { - for (int i = 0; i < corresp.size(); i++) { - corresp.at(i) -= 1; - } - } + if (not quotient_q) + for (int i = 0; i < corresp.size(); i++) corresp.at(i) -= 1; // producing list of simplices std::vector simp_list; @@ -193,7 +200,7 @@ simplicial_complex simplicial_complex::quotient(int char_fun(point_t)) { simp_list.push_back(s_q); } - auto quotient_sc = simplicial_complex(points, simp_list); + simplicial_complex quotient_sc(points, simp_list); quotient_sc.p_impl->quotient_q = quotient_q; return quotient_sc; diff --git a/lib/scomplex/simplicial_complex.hpp b/lib/scomplex/simplicial_complex.hpp index 14c20a9..b1df8fa 100644 --- a/lib/scomplex/simplicial_complex.hpp +++ b/lib/scomplex/simplicial_complex.hpp @@ -10,7 +10,6 @@ using namespace gsimp; class NoChain {}; class simplicial_complex { - private: // implementation details struct impl; std::shared_ptr p_impl; @@ -18,10 +17,11 @@ class simplicial_complex { public: // constructor (no default) simplicial_complex(std::vector&, std::vector&); - + simplicial_complex(const simplicial_complex&); + simplicial_complex& operator=(const simplicial_complex&); // destructor ~simplicial_complex(); - + // std::vector get_points(); std::vector get_level(int); int dimension(); @@ -29,7 +29,7 @@ class simplicial_complex { matrix_t get_boundary_matrix(int); bool is_quotient(); size_t get_index_of_simplex(cell_t); - + // simplicial_complex quotient(int char_fun(point_t)); }; // class simplicial_complex }; // namespace gsimp diff --git a/lib/scomplex/types.hpp b/lib/scomplex/types.hpp index 3cf2783..dd28845 100644 --- a/lib/scomplex/types.hpp +++ b/lib/scomplex/types.hpp @@ -1,4 +1,5 @@ #include +#include #include namespace gsimp { @@ -12,4 +13,7 @@ typedef typename Eigen::SparseVector vector_t; // chains typedef typename std::pair chain_t; + +// characteristic functions +// typedef std::function char_fun_t; }; diff --git a/src/impl_test.cpp b/src/impl_test.cpp index a7403d2..fae9ebc 100644 --- a/src/impl_test.cpp +++ b/src/impl_test.cpp @@ -1,19 +1,28 @@ #include +#include #include #include -int main() { +int my_char(gsimp::point_t point) { return 1; } + +int main(int argc, char* argv[]) { std::vector points_v; std::vector cells_v; std::string filename; - std::cout << "insert file name: "; - std::cin >> filename; + if (argc == 1) { + std::cout << "insert file name: "; + std::cin >> filename; + } else + filename = argv[1]; std::tie(points_v, cells_v) = parse_qhull_file(filename); gsimp::simplicial_complex s_comp(points_v, cells_v); - std::cout << s_comp.dimension(); + + gsimp::simplicial_complex q_comp = s_comp.quotient(my_char); + + std::cout << s_comp.dimension() << "\n"; return 0; } From ff9b5ebc70fdbdb4f7d27977537e61ff9722b3b5 Mon Sep 17 00:00:00 2001 From: crvs Date: Wed, 15 Feb 2017 19:47:34 +0100 Subject: [PATCH 09/24] started implementing the coefficient flow --- lib/scomplex/base_utils.hpp | 6 --- lib/scomplex/coeff_flow.cpp | 50 +++++++++++++++++++ lib/scomplex/qhull_parsing.hpp | 10 ++-- lib/scomplex/simplicial_complex.cpp | 77 +++++++++++++++++++++++++++-- lib/scomplex/simplicial_complex.hpp | 13 ++++- src/impl_test.cpp | 23 ++++++++- 6 files changed, 162 insertions(+), 17 deletions(-) create mode 100644 lib/scomplex/coeff_flow.cpp diff --git a/lib/scomplex/base_utils.hpp b/lib/scomplex/base_utils.hpp index e5471fe..3857270 100644 --- a/lib/scomplex/base_utils.hpp +++ b/lib/scomplex/base_utils.hpp @@ -29,9 +29,3 @@ std::vector dedupe_vec(std::vector& vec) { } return no_reps_list; } - -template -std::unique_ptr new_unique(T* t) { - T* tt = new T(t); - return std::unique_ptr(tt); -} diff --git a/lib/scomplex/coeff_flow.cpp b/lib/scomplex/coeff_flow.cpp new file mode 100644 index 0000000..833d5cd --- /dev/null +++ b/lib/scomplex/coeff_flow.cpp @@ -0,0 +1,50 @@ +#include +#include +#include + +namespace gsimp { + +typedef std::tuple q_elem_t; +typedef std::queue queue_t; + +class out_of_context : std::exception {}; +class no_bounding_chain : std::exception {}; + +int chain_dim(chain_t p) { return std::get<0>(p); } +vector_t chain_rep(chain_t p) { return std::get<1>(p); } + +chain_t coeff_flow(simplicial_complex& s_comp, // + chain_t p, // + simplex_t sigma_0, // + double c_0) { // + + if (chain_dim(p) != s_comp.dim() - 1) throw out_of_context(); + + std::vector seen_sigma(s_comp.level_size(s_comp.dim()), false); + std::vector seen_tau(s_comp.level_size(s_comp.dim() - 1), false); + vector_t b_chain; + queue_t queue; + + cell_t sigma, tau; + double c; + size_t sigma_i, tau_i; + cell_t while (not queue.empty()) { + std::tie(sigma, tau, c) = queue.pop(); + sigma_i = s_comp.get_index_of_simplex(sigma); + tau_i = s_comp.get_index_of_simplex(tau); + if (seen_sigma.at(sigma_i)) { + if (b_chain.coeffRef(sigma_i) != c) throw no_bounding_chain; + } else { + seen_tau.at(tau_i) = true; + seen_sigma.at(sigma_i) = true; + b_chain.coeffRef(sigma_i) = c; + std::vector tau_cofaces; + for (auto coface : s_comp.cofaces(tau)) { + size_t coface_i = s_comp.get_index_of_simplex(coface); + if (coface_i != sigma_i) tau_cofaces.push_back(coface); + } + // need cofaces in index form + } + } +} +}; // namespace gsimp diff --git a/lib/scomplex/qhull_parsing.hpp b/lib/scomplex/qhull_parsing.hpp index d9e5b13..744f9ab 100644 --- a/lib/scomplex/qhull_parsing.hpp +++ b/lib/scomplex/qhull_parsing.hpp @@ -22,10 +22,12 @@ std::vector tokenize(std::string str) { return tokens; } -template -void output_vector(std::vector vec) { - std::copy(vec.begin(), vec.end(), - std::ostream_iterator(std::cout, " ")); +template +void cout_iterator(Iterator first, Iterator last) { + std::copy(first, last, + std::ostream_iterator< + typename std::iterator_traits::value_type>( + std::cout, " ")); std::cout << '\n'; } diff --git a/lib/scomplex/simplicial_complex.cpp b/lib/scomplex/simplicial_complex.cpp index 98813f5..ae93de2 100644 --- a/lib/scomplex/simplicial_complex.cpp +++ b/lib/scomplex/simplicial_complex.cpp @@ -2,10 +2,13 @@ #include #include +#include // for debuging purposes + #include -#include #include -#include + +#include +#include // smart pointers #include #include @@ -50,7 +53,7 @@ struct simplicial_complex::impl { for (auto s : simplex_range) { int d = simplices.dimension(s); simplices.assign_key(s, count.at(d)++); - levels.at(d)->push_back(&s); + levels.at(d)->push_back(new simp_handle(s)); } } @@ -111,8 +114,51 @@ struct simplicial_complex::impl { } return level_cells; } + + simp_handle cell_to_handle(cell_t tau) { + auto sh = simplices.find(tau); + return sh; + } + + cell_t handle_to_cell(simp_handle tau) { + cell_t cell; + for (auto v : simplices.simplex_vertex_range(tau)) cell.push_back(v); + return cell; + } }; // struct impl +int simplicial_complex::boundary_inclusion_index(cell_t c1, cell_t c2) { + return p_impl->boundary_index(p_impl->cell_to_handle(c1), // + p_impl->cell_to_handle(c2)); // +}; +int simplicial_complex::boundary_inclusion_index(int d1, size_t s1, // + int d2, size_t s2) { // + cell_t c1 = index_to_cell(d1, s1); + cell_t c2 = index_to_cell(d2, s2); + return p_impl->boundary_index(p_impl->cell_to_handle(c1), + p_impl->cell_to_handle(c2)); +}; + +std::vector> simplicial_complex::get_cof_and_ind( + cell_t cell) { + std::vector> c_cofaces; + for (auto face : get_cofaces(cell)) + c_cofaces.push_back( // + std::make_pair( // + boundary_inclusion_index(cell, face), face)); // + return c_cofaces; +} + +std::vector> simplicial_complex::get_cof_and_ind_index( + int d, size_t c) { + std::vector> c_cofaces; + for (auto face : get_cofaces_index(d, c)) + c_cofaces.push_back( // + std::make_pair( // + boundary_inclusion_index(d, c, d + 1, face), face)); // + return c_cofaces; +} + int simplicial_complex::get_level_size(int level) { return p_impl->get_level_size(level); } @@ -153,13 +199,36 @@ matrix_t simplicial_complex::get_boundary_matrix(int d) { int simplicial_complex::dimension() { return p_impl->simplices.dimension(); } -size_t simplicial_complex::get_index_of_simplex(cell_t simp) { +cell_t simplicial_complex::index_to_cell(int d, size_t ind) { + auto sh = p_impl->levels.at(d)->at(ind); + return p_impl->handle_to_cell(*sh); +} +size_t simplicial_complex::cell_to_index(cell_t simp) { auto sh = p_impl->simplices.find(simp); return p_impl->simplices.key(sh); } bool simplicial_complex::is_quotient() { return p_impl->quotient_q; } +std::vector simplicial_complex::get_cofaces_index(int d, size_t face) { + std::vector s_cofaces; + impl::simp_handle face_h = *(p_impl->levels.at(d)->at(face)); + cell_t my_cell = p_impl->handle_to_cell(face_h); + // codimension 1 faces + auto range = p_impl->simplices.cofaces_simplex_range(face_h, 1); + for (auto tau : range) s_cofaces.push_back(p_impl->simplices.key(tau)); + return s_cofaces; +} + +std::vector simplicial_complex::get_cofaces(cell_t face) { + std::vector s_cofaces; + auto face_h = p_impl->cell_to_handle(face); + // codimension 1 faces + auto range = p_impl->simplices.cofaces_simplex_range(face_h, 1); + for (auto tau : range) s_cofaces.push_back(p_impl->handle_to_cell(tau)); + return s_cofaces; +} + // TODO: get top dimensional cells from the simplicial complex to add to the // simplex tree (instead of all the simplices). simplicial_complex simplicial_complex::quotient(int char_fun(point_t)) { diff --git a/lib/scomplex/simplicial_complex.hpp b/lib/scomplex/simplicial_complex.hpp index b1df8fa..f2de84e 100644 --- a/lib/scomplex/simplicial_complex.hpp +++ b/lib/scomplex/simplicial_complex.hpp @@ -28,7 +28,18 @@ class simplicial_complex { int get_level_size(int); matrix_t get_boundary_matrix(int); bool is_quotient(); - size_t get_index_of_simplex(cell_t); + // inclusion index + int boundary_inclusion_index(cell_t, cell_t); + int boundary_inclusion_index(int, size_t, int, size_t); + // treating cofaces + std::vector get_cofaces(cell_t); + std::vector get_cofaces_index(int, size_t); + // cofaces and indices jointly + std::vector> get_cof_and_ind(cell_t); + std::vector> get_cof_and_ind_index(int, size_t); + // cells and indices back and forth + cell_t index_to_cell(int, size_t); + size_t cell_to_index(cell_t); // simplicial_complex quotient(int char_fun(point_t)); }; // class simplicial_complex diff --git a/src/impl_test.cpp b/src/impl_test.cpp index fae9ebc..39a4be5 100644 --- a/src/impl_test.cpp +++ b/src/impl_test.cpp @@ -1,5 +1,7 @@ #include +#include +#include #include #include #include @@ -22,7 +24,24 @@ int main(int argc, char* argv[]) { gsimp::simplicial_complex q_comp = s_comp.quotient(my_char); - std::cout << s_comp.dimension() << "\n"; - + std::cout << "complex dimension: " << s_comp.dimension() << "\n"; + for (int i = 0; i <= s_comp.dimension(); ++i) { + std::cout << "level " << i << " has " << s_comp.get_level_size(i) + << " cells\n"; + } + + size_t i; + std::cout << "get 1-cell number: "; + while (std::cin >> i) { + std::vector> my_cofaces = + s_comp.get_cof_and_ind_index(1, i); + for (auto s : my_cofaces) { + std::cout << std::get<0>(s) << ", "; + std::cout << std::get<1>(s) << " ,, "; + } + + std::cout << '\n'; + std::cout << "get 1-cell number: "; + } return 0; } From 3a1f04484a5d986bdc4acecfeb09a0c816c74e9e Mon Sep 17 00:00:00 2001 From: crvs Date: Wed, 15 Feb 2017 21:04:02 +0100 Subject: [PATCH 10/24] continuing coeff flow added facilities for easy coface queries --- lib/scomplex/coeff_flow.cpp | 31 +++++++++++++++++------ lib/scomplex/simplicial_complex.cpp | 39 +++++++++++++++++++++++++++++ lib/scomplex/simplicial_complex.hpp | 13 +++++++--- src/impl_test.cpp | 8 +++--- 4 files changed, 77 insertions(+), 14 deletions(-) diff --git a/lib/scomplex/coeff_flow.cpp b/lib/scomplex/coeff_flow.cpp index 833d5cd..591b4fe 100644 --- a/lib/scomplex/coeff_flow.cpp +++ b/lib/scomplex/coeff_flow.cpp @@ -19,31 +19,46 @@ chain_t coeff_flow(simplicial_complex& s_comp, // double c_0) { // if (chain_dim(p) != s_comp.dim() - 1) throw out_of_context(); - + // 01 + vector_t c_chain(s_comp.level_size(s_comp.dim())); std::vector seen_sigma(s_comp.level_size(s_comp.dim()), false); std::vector seen_tau(s_comp.level_size(s_comp.dim() - 1), false); - vector_t b_chain; + + // 04 + size_t tau_i, sigma_i(s_comp.cell_to_index(sigma_0)); + c_chain.coeffRef(sigma_i) = c_0; + seen_sigma(sigma_i) = true; + + // 05 queue_t queue; + for (cell_t tau : s_comp.cell_boundary(sigma_0)) + queue.push({sigma_0, tau, c_0}); - cell_t sigma, tau; - double c; - size_t sigma_i, tau_i; + // 08 cell_t while (not queue.empty()) { + cell_t sigma, tau; + double c; std::tie(sigma, tau, c) = queue.pop(); - sigma_i = s_comp.get_index_of_simplex(sigma); - tau_i = s_comp.get_index_of_simplex(tau); + sigma_i = s_comp.cell_to_index(sigma); + tau_i = s_comp.cell_to_index(tau); + + // 10 if (seen_sigma.at(sigma_i)) { if (b_chain.coeffRef(sigma_i) != c) throw no_bounding_chain; } else { + // 14 seen_tau.at(tau_i) = true; seen_sigma.at(sigma_i) = true; b_chain.coeffRef(sigma_i) = c; + + // 15 + // + // std::vector tau_cofaces; for (auto coface : s_comp.cofaces(tau)) { size_t coface_i = s_comp.get_index_of_simplex(coface); if (coface_i != sigma_i) tau_cofaces.push_back(coface); } - // need cofaces in index form } } } diff --git a/lib/scomplex/simplicial_complex.cpp b/lib/scomplex/simplicial_complex.cpp index ae93de2..7846df3 100644 --- a/lib/scomplex/simplicial_complex.cpp +++ b/lib/scomplex/simplicial_complex.cpp @@ -115,6 +115,12 @@ struct simplicial_complex::impl { return level_cells; } + simp_handle index_to_handle(int d, size_t tau) { + return *(levels.at(d)->at(tau)); + } + + size_t handle_to_index(simp_handle tau) { return simplices.key(tau); } + simp_handle cell_to_handle(cell_t tau) { auto sh = simplices.find(tau); return sh; @@ -127,6 +133,39 @@ struct simplicial_complex::impl { } }; // struct impl +std::vector simplicial_complex::cell_boundary_index(cell_t cell) { + impl::simp_handle simp = p_impl->cell_to_handle(cell); + auto c_boundary = p_impl->simplices.boundary_simplex_range(simp); + std::vector s_boundary; + for (auto c : c_boundary) s_boundary.push_back(p_impl->handle_to_index(c)); + return s_boundary; +} + +std::vector simplicial_complex::cell_boundary_index(int d, + size_t cell) { + impl::simp_handle simp = p_impl->index_to_handle(d, cell); + auto c_boundary = p_impl->simplices.boundary_simplex_range(simp); + std::vector s_boundary; + for (auto c : c_boundary) s_boundary.push_back(p_impl->handle_to_index(c)); + return s_boundary; +} + +std::vector simplicial_complex::cell_boundary(cell_t cell) { + impl::simp_handle simp = p_impl->cell_to_handle(cell); + auto c_boundary = p_impl->simplices.boundary_simplex_range(simp); + std::vector s_boundary; + for (auto c : c_boundary) s_boundary.push_back(p_impl->handle_to_cell(c)); + return s_boundary; +} + +std::vector simplicial_complex::cell_boundary(int d, size_t cell) { + impl::simp_handle simp = p_impl->index_to_handle(d, cell); + auto c_boundary = p_impl->simplices.boundary_simplex_range(simp); + std::vector s_boundary; + for (auto c : c_boundary) s_boundary.push_back(p_impl->handle_to_cell(c)); + return s_boundary; +} + int simplicial_complex::boundary_inclusion_index(cell_t c1, cell_t c2) { return p_impl->boundary_index(p_impl->cell_to_handle(c1), // p_impl->cell_to_handle(c2)); // diff --git a/lib/scomplex/simplicial_complex.hpp b/lib/scomplex/simplicial_complex.hpp index f2de84e..1d8edd7 100644 --- a/lib/scomplex/simplicial_complex.hpp +++ b/lib/scomplex/simplicial_complex.hpp @@ -21,13 +21,20 @@ class simplicial_complex { simplicial_complex& operator=(const simplicial_complex&); // destructor ~simplicial_complex(); - // + // basic info std::vector get_points(); - std::vector get_level(int); int dimension(); + bool is_quotient(); + // level-wise info + std::vector get_level(int); int get_level_size(int); + // cell boundaries + std::vector cell_boundary_index(cell_t); + std::vector cell_boundary_index(int, size_t); + std::vector cell_boundary(cell_t); + std::vector cell_boundary(int, size_t); + // boundary matrices matrix_t get_boundary_matrix(int); - bool is_quotient(); // inclusion index int boundary_inclusion_index(cell_t, cell_t); int boundary_inclusion_index(int, size_t, int, size_t); diff --git a/src/impl_test.cpp b/src/impl_test.cpp index 39a4be5..8e7b667 100644 --- a/src/impl_test.cpp +++ b/src/impl_test.cpp @@ -34,14 +34,16 @@ int main(int argc, char* argv[]) { std::cout << "get 1-cell number: "; while (std::cin >> i) { std::vector> my_cofaces = - s_comp.get_cof_and_ind_index(1, i); + s_comp.get_cof_and_ind_index(2, i); for (auto s : my_cofaces) { std::cout << std::get<0>(s) << ", "; std::cout << std::get<1>(s) << " ,, "; } - std::cout << '\n'; - std::cout << "get 1-cell number: "; + std::vector faces = s_comp.cell_boundary_index(2, i); + std::copy(faces.begin(), faces.end(), + std::ostream_iterator(std::cout, " ")); + std::cout << "\nget 1-cell number: "; } return 0; } From 610616047dede9e5b6134d4587f20a9476cc5713 Mon Sep 17 00:00:00 2001 From: crvs Date: Thu, 16 Feb 2017 10:34:41 +0100 Subject: [PATCH 11/24] switched to clang --- CMakeLists.txt | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index e54d9c3..4b9d27c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -21,6 +21,26 @@ include(${CGAL_USE_FILE}) # end_requirements +# Using CLANG instead of gcc for reasons +set(CMAKE_C_COMPILER "/usr/bin/clang") +set(CMAKE_C_FLAGS "-Wall -std=c99") +set(CMAKE_C_FLAGS_DEBUG "-g") +set(CMAKE_C_FLAGS_MINSIZEREL "-Os -DNDEBUG") +set(CMAKE_C_FLAGS_RELEASE "-O4 -DNDEBUG") +set(CMAKE_C_FLAGS_RELWITHDEBINFO "-O2 -g") + +set(CMAKE_CXX_COMPILER "/usr/bin/clang++") +set(CMAKE_CXX_FLAGS "-Wall") +set(CMAKE_CXX_FLAGS_DEBUG "-g") +set(CMAKE_CXX_FLAGS_MINSIZEREL "-Os -DNDEBUG") +set(CMAKE_CXX_FLAGS_RELEASE "-O4 -DNDEBUG") +set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O2 -g") + +set(CMAKE_AR "/usr/bin/llvm-ar") +set(CMAKE_LINKER "/usr/bin/llvm-ld") +set(CMAKE_NM "/usr/bin/llvm-nm") +set(CMAKE_OBJDUMP "/usr/bin/llvm-objdump") +set(CMAKE_RANLIB "/usr/bin/llvm-ranlib") set(CMAKE_BUILD_TYPE Debug) include_directories("./lib") # where the header files are From 3d6f58f3ac0ff14014390a6ae5345810664f1fec Mon Sep 17 00:00:00 2001 From: crvs Date: Thu, 16 Feb 2017 17:34:39 +0100 Subject: [PATCH 12/24] started refactoring of path snapper --- lib/scomplex/path-snapping.hpp | 28 ++++++++++ lib/scomplex/path_snapper.cpp | 78 ++++++++++++++++++++++++++++ lib/scomplex/path_snapper.hpp | 93 +++++++++------------------------- lib/scomplex/types.hpp | 4 +- 4 files changed, 132 insertions(+), 71 deletions(-) create mode 100644 lib/scomplex/path-snapping.hpp create mode 100644 lib/scomplex/path_snapper.cpp diff --git a/lib/scomplex/path-snapping.hpp b/lib/scomplex/path-snapping.hpp new file mode 100644 index 0000000..ca5bce1 --- /dev/null +++ b/lib/scomplex/path-snapping.hpp @@ -0,0 +1,28 @@ +#pragma once + +namespace snap { + +#include +#include +#include +#include +#include + +typedef Eigen::SparseVector vector_t; + +class path_snapper { + private: + struct impl; + std::shared_ptr p_impl; + + public: + path_snapper(std::vector& pts, std::vector& cells); + path_snapper(simplicial::simplicial_complex&); + ~path_snapper(); + path_snapper(path_snapper& other); + path_snapper& operator=(path_snapper& other); + + std::vector snap_path(std::vector path); + vector_t get_chain_vector(std::vector path); +}; +}; diff --git a/lib/scomplex/path_snapper.cpp b/lib/scomplex/path_snapper.cpp new file mode 100644 index 0000000..dbcf900 --- /dev/null +++ b/lib/scomplex/path_snapper.cpp @@ -0,0 +1,78 @@ +#pragma once + +#include +#include + +namespace gsimp { + +#include +#include +#include + +struct path_snapper::impl { + tree_t point_tree; // defined in nn_utils.hpp + Graph vertex_graph; // defined in graph_utils.hpp + simplicial::SimplicialComplex* s_comp; + bool owner; + + typedef std::vector, int>> geo_chain_t; + geo_chain_t get_chain_rep(std::vector path) { + auto vertex_path = snap_path(path); + geo_chain_t chain_rep; + for ( // + auto it = vertex_path.begin(); // + std::next(it) != vertex_path.end(); // + ++it) { + size_t s(*it), t(*(std::next(it))); + if (s < t) { + chain_rep.push_back( // + std::make_pair, int>( + std::vector{s, t}, 1)); + } else { + chain_rep.push_back( // + std::make_pair, int>( + std::vector{s, t}, -1)); + } + } + return chain_rep; + } +} + +public : path_snapper(simplicial::SimplicialComplex& sc) { + owner = false; + s_comp = ≻ + make_tree(point_tree, s_comp->points); + vertex_graph = calculate_one_skelleton_graph(*s_comp); +} + +path_snapper(std::vector& pts, std::vector& cells) { + owner = true; + s_comp = new simplicial::SimplicialComplex(pts, cells); + make_tree(point_tree, s_comp->points); + vertex_graph = calculate_one_skelleton_graph(*s_comp); +} + +~path_snapper() { + if (owner) delete s_comp; +} + +std::vector snap_path(std::vector path) { + std::vector way_points = snap_points_to_indexes(point_tree, path); + std::vector snapped_path = complete_path(vertex_graph, way_points); + return snapped_path; +} + +vector_t get_chain_vector(std::vector path) { + auto chain_rep = get_chain_rep(path); + vector_t vector_rep(s_comp->get_dimension(1)); + for (auto chain_el : chain_rep) { + simplex_t simp = std::get<0>(chain_el); + int val = std::get<1>(chain_el); + size_t index = s_comp->get_simplex_index(simp); + vector_rep.coeffRef(index) = val; + } + return vector_rep; +} +}; +} +; diff --git a/lib/scomplex/path_snapper.hpp b/lib/scomplex/path_snapper.hpp index ca1d9b2..0c8cee0 100644 --- a/lib/scomplex/path_snapper.hpp +++ b/lib/scomplex/path_snapper.hpp @@ -1,80 +1,35 @@ #pragma once -namespace snap { +#include +#include -#include -#include -#include -#include +namespace gsimp { -typedef Eigen::SparseVector vector_t; +#include +#include +#include class path_snapper { private: - tree_t point_tree; // defined in nn_utils.hpp - Graph vertex_graph; // defined in graph_utils.hpp - simplicial::SimplicialComplex* s_comp; - bool owner; - - typedef std::vector, int>> geo_chain_t; - geo_chain_t get_chain_rep(std::vector path) { - auto vertex_path = snap_path(path); - geo_chain_t chain_rep; - for ( // - auto it = vertex_path.begin(); // - std::next(it) != vertex_path.end(); // - ++it) { - size_t s(*it), t(*(std::next(it))); - if (s < t) { - chain_rep.push_back( // - std::make_pair, int>( - std::vector{s, t}, 1)); - } else { - chain_rep.push_back( // - std::make_pair, int>( - std::vector{s, t}, -1)); - } - } - return chain_rep; - } + struct impl; + std::shared_ptr p_impl; public: - path_snapper(simplicial::SimplicialComplex& sc) { - owner = false; - s_comp = ≻ - make_tree(point_tree, s_comp->points); - vertex_graph = calculate_one_skelleton_graph(*s_comp); - } - - path_snapper(std::vector& pts, std::vector& cells) { - owner = true; - s_comp = new simplicial::SimplicialComplex(pts, cells); - make_tree(point_tree, s_comp->points); - vertex_graph = calculate_one_skelleton_graph(*s_comp); - } - - ~path_snapper() { - if (owner) delete s_comp; - } - - std::vector snap_path(std::vector path) { - std::vector way_points = - snap_points_to_indexes(point_tree, path); - std::vector snapped_path = - complete_path(vertex_graph, way_points); - return snapped_path; - } - - vector_t get_chain_vector(std::vector path) { - auto chain_rep = get_chain_rep(path); - vector_t vector_rep(s_comp->get_dimension(1)); - for (auto chain_el : chain_rep) { - simplex_t simp = std::get<0>(chain_el); - int val = std::get<1>(chain_el); - size_t index = s_comp->get_simplex_index(simp); - vector_rep.coeffRef(index) = val; - } - return vector_rep; - } + // constructors and such + path_snapper(std::vector&, std::vector&); + path_snapper(simplicial::simplicial_complex&); + ~path_snapper(); + path_snapper(path_snapper&); + path_snapper& operator=(path_snapper&); + // the things we want to do + std::vector snap_path_to_points(std::vector); + std::vector snap_path_to_indices(std::vector); + vector_t get_chain_vector(std::vector); + chain_t snap_path_to_chain(std::vector); + // interconversion + std::vector index_sequence_to_point(std::vector); + std::vector point_sequence_to_index(std::vector); + chain_t index_sequence_to_chain(std::vector); + chain_t point_sequence_to_chain(std::vector); }; }; diff --git a/lib/scomplex/types.hpp b/lib/scomplex/types.hpp index dd28845..511b445 100644 --- a/lib/scomplex/types.hpp +++ b/lib/scomplex/types.hpp @@ -1,6 +1,6 @@ -#include -#include #include +#include +#include namespace gsimp { // geometric types From 3b5994a65c65c5d46f0f0c33ad137f6d2a48222f Mon Sep 17 00:00:00 2001 From: crvs Date: Thu, 16 Feb 2017 18:57:08 +0100 Subject: [PATCH 13/24] finished coeff flow implementation --- CMakeLists.txt | 12 +++-- lib/scomplex/coeff_flow.cpp | 72 ++++++++++++++++++++--------- lib/scomplex/qhull_parsing.hpp | 6 +-- lib/scomplex/simplicial_complex.cpp | 39 ++++++++++------ lib/scomplex/simplicial_complex.hpp | 17 +++---- src/impl_test.cpp | 3 ++ 6 files changed, 92 insertions(+), 57 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 4b9d27c..702c65d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -30,11 +30,13 @@ set(CMAKE_C_FLAGS_RELEASE "-O4 -DNDEBUG") set(CMAKE_C_FLAGS_RELWITHDEBINFO "-O2 -g") set(CMAKE_CXX_COMPILER "/usr/bin/clang++") -set(CMAKE_CXX_FLAGS "-Wall") +set(CMAKE_CXX_FLAGS "-Wall -std=c++11") set(CMAKE_CXX_FLAGS_DEBUG "-g") set(CMAKE_CXX_FLAGS_MINSIZEREL "-Os -DNDEBUG") set(CMAKE_CXX_FLAGS_RELEASE "-O4 -DNDEBUG") set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O2 -g") +# Linker flags for +# set(CMAKE_CXX_FLAGS "-lspatialindex_c -lspatialindex ${CMAKE_CXX_FLAGS}") set(CMAKE_AR "/usr/bin/llvm-ar") set(CMAKE_LINKER "/usr/bin/llvm-ld") @@ -43,14 +45,14 @@ set(CMAKE_OBJDUMP "/usr/bin/llvm-objdump") set(CMAKE_RANLIB "/usr/bin/llvm-ranlib") set(CMAKE_BUILD_TYPE Debug) + include_directories("./lib") # where the header files are -set(CMAKE_CXX_FLAGS "-std=c++11 ${CMAKE_CXX_FLAGS}") -set(CMAKE_CXX_FLAGS "-lspatialindex_c -lspatialindex ${CMAKE_CXX_FLAGS}") -add_library(test SHARED "./lib/scomplex/simplicial_complex.cpp") +add_library(scomplex SHARED "./lib/scomplex/simplicial_complex.cpp") +add_library(coeffflow SHARED "./lib/scomplex/coeff_flow.cpp") add_executable(impl_test "src/impl_test.cpp") -target_link_libraries(impl_test test) +target_link_libraries(impl_test scomplex) # set(STUPID_SOURCES "./src/main.cpp") # add_executable(my_stupid_test ${STUPID_SOURCES}) diff --git a/lib/scomplex/coeff_flow.cpp b/lib/scomplex/coeff_flow.cpp index 591b4fe..29c825a 100644 --- a/lib/scomplex/coeff_flow.cpp +++ b/lib/scomplex/coeff_flow.cpp @@ -1,65 +1,93 @@ #include #include +#include +#include #include namespace gsimp { typedef std::tuple q_elem_t; -typedef std::queue queue_t; +typedef std::queue queue_t; class out_of_context : std::exception {}; class no_bounding_chain : std::exception {}; -int chain_dim(chain_t p) { return std::get<0>(p); } -vector_t chain_rep(chain_t p) { return std::get<1>(p); } +int& chain_dim(chain_t& p) { return std::get<0>(p); } +vector_t& chain_rep(chain_t& p) { return std::get<1>(p); } +double& chain_val(chain_t& p, size_t i) { return std::get<1>(p).coeffRef(i); } chain_t coeff_flow(simplicial_complex& s_comp, // chain_t p, // - simplex_t sigma_0, // + cell_t sigma_0, // double c_0) { // - if (chain_dim(p) != s_comp.dim() - 1) throw out_of_context(); + if (chain_dim(p) != s_comp.dimension() - 1) throw out_of_context(); // 01 - vector_t c_chain(s_comp.level_size(s_comp.dim())); - std::vector seen_sigma(s_comp.level_size(s_comp.dim()), false); - std::vector seen_tau(s_comp.level_size(s_comp.dim() - 1), false); + chain_t c_chain(s_comp.dimension(), + vector_t(s_comp.get_level_size(s_comp.dimension()))); + std::vector seen_sigma(s_comp.get_level_size(s_comp.dimension()), + false); + std::vector seen_tau(s_comp.get_level_size(s_comp.dimension() - 1), + false); // 04 size_t tau_i, sigma_i(s_comp.cell_to_index(sigma_0)); - c_chain.coeffRef(sigma_i) = c_0; - seen_sigma(sigma_i) = true; + chain_val(c_chain, sigma_i) = c_0; + seen_sigma.at(sigma_i) = true; // 05 queue_t queue; for (cell_t tau : s_comp.cell_boundary(sigma_0)) - queue.push({sigma_0, tau, c_0}); + queue.emplace(sigma_0, tau, c_0); // 08 - cell_t while (not queue.empty()) { + while (not queue.empty()) { cell_t sigma, tau; double c; - std::tie(sigma, tau, c) = queue.pop(); + std::tie(sigma, tau, c) = queue.front(); + queue.pop(); sigma_i = s_comp.cell_to_index(sigma); tau_i = s_comp.cell_to_index(tau); // 10 - if (seen_sigma.at(sigma_i)) { - if (b_chain.coeffRef(sigma_i) != c) throw no_bounding_chain; - } else { + if (seen_sigma.at(sigma_i) && chain_val(c_chain, sigma_i) != c) + throw no_bounding_chain(); + else if (not seen_sigma.at(sigma_i)) { // 14 seen_tau.at(tau_i) = true; seen_sigma.at(sigma_i) = true; - b_chain.coeffRef(sigma_i) = c; + chain_val(c_chain, sigma_i) = c; // 15 - // - // - std::vector tau_cofaces; - for (auto coface : s_comp.cofaces(tau)) { - size_t coface_i = s_comp.get_index_of_simplex(coface); + std::vector> tau_cofaces; + for (auto coface : s_comp.get_cof_and_ind(tau)) { + size_t coface_i = s_comp.cell_to_index(std::get<1>(coface)); if (coface_i != sigma_i) tau_cofaces.push_back(coface); } + + // 16 + double predicted_bdry = + s_comp.boundary_inclusion_index(tau, sigma) * c; + if (tau_cofaces.size() == 0 && // + predicted_bdry != chain_val(p, tau_i)) // + throw no_bounding_chain(); // 18 + else if (tau_cofaces.size() != 0) { + int sigma_p_ind; + cell_t sigma_p; + std::tie(sigma_p_ind, sigma_p) = tau_cofaces[0]; // 20 + double c_p = + (sigma_p_ind * + (chain_val(p, tau_i) - + s_comp.boundary_inclusion_index(tau, sigma) * c)); + for (cell_t tau_p : s_comp.cell_boundary(sigma_p)) { + if (not seen_tau.at(s_comp.cell_to_index(tau_p))) { + queue.emplace(sigma_p, tau_p, c_p); + seen_tau.at(tau_i) = true; + } + } + } } } + return c_chain; } }; // namespace gsimp diff --git a/lib/scomplex/qhull_parsing.hpp b/lib/scomplex/qhull_parsing.hpp index 744f9ab..446b2db 100644 --- a/lib/scomplex/qhull_parsing.hpp +++ b/lib/scomplex/qhull_parsing.hpp @@ -41,11 +41,7 @@ std::pair, std::vector> parse_qhull_file( std::string line; - int current_line_number{0}, dimensionality; - int num_points, num_cells; - int cur_point_num, cur_cell_num; - - int line_number{0}, total_points{0}, total_cells{0}; + int dimensionality, total_points{0}, total_cells{0}; std::vector points; std::vector cells; diff --git a/lib/scomplex/simplicial_complex.cpp b/lib/scomplex/simplicial_complex.cpp index 7846df3..f57b4b7 100644 --- a/lib/scomplex/simplicial_complex.cpp +++ b/lib/scomplex/simplicial_complex.cpp @@ -34,7 +34,7 @@ struct simplicial_complex::impl { impl() : quotient_q(false){}; impl(std::vector& arg_points, std::vector& arg_tris) - : quotient_q(false), points(arg_points) { + : points(arg_points), quotient_q(false) { // create the simplex tree for (auto tri : arg_tris) { simplices.insert_simplex_and_subfaces(dedupe_vec(tri)); @@ -70,7 +70,6 @@ struct simplicial_complex::impl { auto end_1 = simplices.simplex_vertex_range(s_1).end(); // auto it_2 = simplices.simplex_vertex_range(s_2).begin(); - auto end_2 = simplices.simplex_vertex_range(s_2).end(); // int orient = 0; int unmatch = 0; @@ -133,13 +132,31 @@ struct simplicial_complex::impl { } }; // struct impl -std::vector simplicial_complex::cell_boundary_index(cell_t cell) { +std::vector> simplicial_complex::get_bdry_and_ind( + cell_t cell) { impl::simp_handle simp = p_impl->cell_to_handle(cell); auto c_boundary = p_impl->simplices.boundary_simplex_range(simp); - std::vector s_boundary; - for (auto c : c_boundary) s_boundary.push_back(p_impl->handle_to_index(c)); - return s_boundary; -} + std::vector> boundary_and_indices; + for (auto face : c_boundary) + boundary_and_indices.push_back( // + std::make_pair( // + p_impl->boundary_index(face, simp), // + p_impl->handle_to_cell(face))); // + return boundary_and_indices; +}; + +std::vector> simplicial_complex::get_bdry_and_ind_index( + int d, size_t cell) { + impl::simp_handle simp = p_impl->index_to_handle(d, cell); + auto c_boundary = p_impl->simplices.boundary_simplex_range(simp); + std::vector> boundary_and_indices; + for (auto face : c_boundary) + boundary_and_indices.push_back( // + std::make_pair( // + p_impl->boundary_index(face, simp), // + p_impl->handle_to_index(face))); // + return boundary_and_indices; +}; std::vector simplicial_complex::cell_boundary_index(int d, size_t cell) { @@ -158,14 +175,6 @@ std::vector simplicial_complex::cell_boundary(cell_t cell) { return s_boundary; } -std::vector simplicial_complex::cell_boundary(int d, size_t cell) { - impl::simp_handle simp = p_impl->index_to_handle(d, cell); - auto c_boundary = p_impl->simplices.boundary_simplex_range(simp); - std::vector s_boundary; - for (auto c : c_boundary) s_boundary.push_back(p_impl->handle_to_cell(c)); - return s_boundary; -} - int simplicial_complex::boundary_inclusion_index(cell_t c1, cell_t c2) { return p_impl->boundary_index(p_impl->cell_to_handle(c1), // p_impl->cell_to_handle(c2)); // diff --git a/lib/scomplex/simplicial_complex.hpp b/lib/scomplex/simplicial_complex.hpp index 1d8edd7..1fa21af 100644 --- a/lib/scomplex/simplicial_complex.hpp +++ b/lib/scomplex/simplicial_complex.hpp @@ -5,8 +5,6 @@ namespace gsimp { -using namespace gsimp; - class NoChain {}; class simplicial_complex { @@ -28,22 +26,21 @@ class simplicial_complex { // level-wise info std::vector get_level(int); int get_level_size(int); - // cell boundaries - std::vector cell_boundary_index(cell_t); - std::vector cell_boundary_index(int, size_t); - std::vector cell_boundary(cell_t); - std::vector cell_boundary(int, size_t); - // boundary matrices - matrix_t get_boundary_matrix(int); // inclusion index int boundary_inclusion_index(cell_t, cell_t); int boundary_inclusion_index(int, size_t, int, size_t); + // cell boundaries + std::vector cell_boundary(cell_t); + std::vector cell_boundary_index(int, size_t); + std::vector> get_bdry_and_ind(cell_t); + std::vector> get_bdry_and_ind_index(int, size_t); // treating cofaces std::vector get_cofaces(cell_t); std::vector get_cofaces_index(int, size_t); - // cofaces and indices jointly std::vector> get_cof_and_ind(cell_t); std::vector> get_cof_and_ind_index(int, size_t); + // boundary matrices + matrix_t get_boundary_matrix(int); // cells and indices back and forth cell_t index_to_cell(int, size_t); size_t cell_to_index(cell_t); diff --git a/src/impl_test.cpp b/src/impl_test.cpp index 8e7b667..6e82b58 100644 --- a/src/impl_test.cpp +++ b/src/impl_test.cpp @@ -8,6 +8,8 @@ int my_char(gsimp::point_t point) { return 1; } +int& do_this(int& a) { return a; } + int main(int argc, char* argv[]) { std::vector points_v; std::vector cells_v; @@ -45,5 +47,6 @@ int main(int argc, char* argv[]) { std::ostream_iterator(std::cout, " ")); std::cout << "\nget 1-cell number: "; } + return 0; } From 4710bd02aad3e5e16ff123fc40525325978537de Mon Sep 17 00:00:00 2001 From: crvs Date: Fri, 17 Feb 2017 10:27:59 +0100 Subject: [PATCH 14/24] notes on drawing --- lib/scomplex/chain-calc.hpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/lib/scomplex/chain-calc.hpp b/lib/scomplex/chain-calc.hpp index 0e4dc47..dd028c4 100644 --- a/lib/scomplex/chain-calc.hpp +++ b/lib/scomplex/chain-calc.hpp @@ -1,12 +1,12 @@ #pragma once +#include +#include +#include +#include #include #include -#include #include -#include -#include -#include #include #include @@ -87,6 +87,10 @@ bool equals(vector_t vec1, vector_t vec2) { return true; } +// TODO: drowing facilities +// - draw chain +// - draw 1-chain with arrows +// - put it into its own header file void bounding_chain::draw_chain(std::string filename, chain_t chain) { typedef boost::geometry::model::d2::point_xy point_type; typedef boost::geometry::model::ring polygon_type; From 1c333a23d003ea30a5296f0ecf8ed7c096d9c0bb Mon Sep 17 00:00:00 2001 From: crvs Date: Fri, 17 Feb 2017 12:13:27 +0100 Subject: [PATCH 15/24] Refactoring - using clang-3.8 - moving chain_t utilities from coeff_flow.hpp to types.hpp - guarded types.hpp - added an empty chain creator to simplicial_complex - refactored graph_utils to comply with the changes made - fixed path_snapper (header) constructor definitions and made the get_chain return a chain - completed path_snapper implementation --- CMakeLists.txt | 5 +- lib/scomplex/base_utils.hpp | 31 ------ lib/scomplex/coeff_flow.cpp | 4 - lib/scomplex/graph_utils.hpp | 102 ++++++++++--------- lib/scomplex/path_snapper.cpp | 147 ++++++++++++++++++---------- lib/scomplex/path_snapper.hpp | 12 +-- lib/scomplex/simplicial_complex.cpp | 20 +++- lib/scomplex/simplicial_complex.hpp | 2 + lib/scomplex/types.hpp | 7 ++ 9 files changed, 178 insertions(+), 152 deletions(-) delete mode 100644 lib/scomplex/base_utils.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 702c65d..f9577bb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,14 +22,14 @@ include(${CGAL_USE_FILE}) # Using CLANG instead of gcc for reasons -set(CMAKE_C_COMPILER "/usr/bin/clang") +set(CMAKE_C_COMPILER "/usr/bin/clang-3.8") set(CMAKE_C_FLAGS "-Wall -std=c99") set(CMAKE_C_FLAGS_DEBUG "-g") set(CMAKE_C_FLAGS_MINSIZEREL "-Os -DNDEBUG") set(CMAKE_C_FLAGS_RELEASE "-O4 -DNDEBUG") set(CMAKE_C_FLAGS_RELWITHDEBINFO "-O2 -g") -set(CMAKE_CXX_COMPILER "/usr/bin/clang++") +set(CMAKE_CXX_COMPILER "/usr/bin/clang++-3.8") set(CMAKE_CXX_FLAGS "-Wall -std=c++11") set(CMAKE_CXX_FLAGS_DEBUG "-g") set(CMAKE_CXX_FLAGS_MINSIZEREL "-Os -DNDEBUG") @@ -50,6 +50,7 @@ include_directories("./lib") # where the header files are add_library(scomplex SHARED "./lib/scomplex/simplicial_complex.cpp") add_library(coeffflow SHARED "./lib/scomplex/coeff_flow.cpp") +add_library(pathsnap SHARED "./lib/scomplex/path_snapper.cpp") add_executable(impl_test "src/impl_test.cpp") target_link_libraries(impl_test scomplex) diff --git a/lib/scomplex/base_utils.hpp b/lib/scomplex/base_utils.hpp deleted file mode 100644 index 3857270..0000000 --- a/lib/scomplex/base_utils.hpp +++ /dev/null @@ -1,31 +0,0 @@ -#pragma once - -#include -#include -#include - -std::list dedupe_list(std::list& list) { - // O(n log(n)) - std::set no_reps; - std::list no_reps_list; - for (size_t el : list) { - no_reps.insert(el); - } - for (size_t el : no_reps) { - no_reps_list.push_back(el); - } - return no_reps_list; -} - -std::vector dedupe_vec(std::vector& vec) { - // O(n log(n)) - std::set no_reps; - std::vector no_reps_list; - for (size_t el : vec) { - no_reps.insert(el); - } - for (size_t el : no_reps) { - no_reps_list.push_back(el); - } - return no_reps_list; -} diff --git a/lib/scomplex/coeff_flow.cpp b/lib/scomplex/coeff_flow.cpp index 29c825a..44f084b 100644 --- a/lib/scomplex/coeff_flow.cpp +++ b/lib/scomplex/coeff_flow.cpp @@ -12,10 +12,6 @@ typedef std::queue queue_t; class out_of_context : std::exception {}; class no_bounding_chain : std::exception {}; -int& chain_dim(chain_t& p) { return std::get<0>(p); } -vector_t& chain_rep(chain_t& p) { return std::get<1>(p); } -double& chain_val(chain_t& p, size_t i) { return std::get<1>(p).coeffRef(i); } - chain_t coeff_flow(simplicial_complex& s_comp, // chain_t p, // cell_t sigma_0, // diff --git a/lib/scomplex/graph_utils.hpp b/lib/scomplex/graph_utils.hpp index e13346c..becf276 100644 --- a/lib/scomplex/graph_utils.hpp +++ b/lib/scomplex/graph_utils.hpp @@ -1,13 +1,8 @@ #pragma once #include -// reverse -#include - -// to compute norms -#include "Eigen/Eigen" - -#include "scomplex/Simplicial_Complex.h" +#include +#include // graph libraries #include "boost/graph/adjacency_list.hpp" @@ -15,63 +10,68 @@ #include "boost/graph/named_function_params.hpp" #include "boost/property_map/vector_property_map.hpp" +namespace gsimp { + +/* +typedefs: + graph_t +functions: + calculate_one_skelleton_graph(simplicial_complex& s_comp ) + complete_path(const graph_t& g, std::vector vec) + shortest_path(const graph_t& g, vertex_t s, vertex_t t) +*/ + typedef typename boost::adjacency_list< // boost::vecS, // boost::vecS, // boost::undirectedS, // boost::no_property, // vertex property boost::property> // - Graph; + graph_t; -typedef typename boost::graph_traits::vertex_descriptor vertex_t; +typedef typename boost::graph_traits::vertex_descriptor vertex_t; -typedef std::vector point_t; -typedef std::vector simplex_t; - -Graph calculate_one_skelleton_graph(simplicial::SimplicialComplex& scomp // - ) { - int num_points{scomp.get_dimension(0)}; - int num_edges{scomp.get_dimension(1)}; +graph_t calculate_one_skelleton_graph(simplicial_complex& s_comp // + ) { + size_t num_points(s_comp.get_level_size(0)); + size_t num_edges(s_comp.get_level_size(1)); std::vector> edges(num_edges); std::vector weights(num_edges); - for (simplex_t edge : scomp.get_level(1)) { + for (cell_t edge : s_comp.get_level(1)) { int src_index = edge.at(0); int trg_index = edge.at(1); - point_t src_point = scomp.points.at(edge.at(0)); - point_t trg_point = scomp.points.at(edge.at(1)); + point_t src_point = s_comp.get_points().at(edge.at(0)); + point_t trg_point = s_comp.get_points().at(edge.at(1)); edges.push_back(std::make_pair(src_index, trg_index)); - Eigen::VectorXd diff(src_point.size()); - - for (int i = 0; i < src_point.size(); i++) - diff(i) = src_point.at(i) - trg_point.at(i); - - weights.push_back(diff.norm()); + double norm = 0; + for (int i = 0; i < src_point.size(); ++i) + norm += pow(src_point.at(i) - trg_point.at(i), 2); + weights.push_back(sqrt(norm)); } - const Graph one_skelleton(edges.begin(), // - edges.end(), // - weights.begin(), // - num_points); // - + const graph_t one_skelleton(edges.begin(), // + edges.end(), // + weights.begin(), // + num_points); // return one_skelleton; } /** - * @brief find the shortest path between two vertices in a graph. - * - * @param g: graph (passed by ref) - * @param s: source vertex (type: decltype(g)::vertex_descriptor) - * @param t: source vertex (type: decltype(g)::vertex_descriptor) - * - * @returns: a vector of vertices - * (type: std::vector) - */ -std::vector shortest_path(const Graph& g, vertex_t s, vertex_t t) { +* @brief find the shortest path between two vertices in a graph. +* +* @param g: graph (passed by ref) +* @param s: source vertex (type: decltype(g)::vertex_descriptor) +* @param t: source vertex (type: decltype(g)::vertex_descriptor) +* +* @returns: a vector of vertices +* (type: std::vector) +*/ +std::vector shortest_path(const graph_t& g, vertex_t s, vertex_t t) { // make a predecessor map boost::vector_property_map predecessors(boost::num_vertices(g)); boost::dijkstra_shortest_paths(g, s, @@ -79,18 +79,15 @@ std::vector shortest_path(const Graph& g, vertex_t s, vertex_t t) { // construct the path (going backwards from the target to the source) vertex_t it = t; std::vector s_t_path{}; - // WARNING: the push_back and reverse strategy may be more efficient - // - // use a do-while block, we want the guard to be checked afterwards do { - // push to the front - s_t_path.insert(s_t_path.begin(), it); - it = predecessors[it]; + s_t_path.push_back(it); } while (it != s); + std::reverse(s_t_path.begin(), s_t_path.end()); return s_t_path; } -std::vector complete_path(const Graph& g, std::vector vec) { +std::vector complete_path(const graph_t& g, + std::vector vec) { std::vector::iterator s{vec.begin()}; std::vector::iterator t = std::next(s); @@ -98,15 +95,14 @@ std::vector complete_path(const Graph& g, std::vector vec) { std::vector full_path; full_path.push_back(*s); while (t != vec.end()) { - // get the next vertex - // calculate and add the path between s and t (contains vertices other - // than s) to the full path + // get the next vertex calculate and add the path between s and t + // (contains vertices other than s) to the full path std::vector path_portion{shortest_path(g, *s, *t)}; - for (vertex_t vert : path_portion) { - full_path.push_back(vert); - } + // quick vector merging + full_path.insert(full_path.end(),path_portion.begin(),path_portion.end()); std::advance(s, 1); std::advance(t, 1); } return full_path; } +}; diff --git a/lib/scomplex/path_snapper.cpp b/lib/scomplex/path_snapper.cpp index dbcf900..1a62637 100644 --- a/lib/scomplex/path_snapper.cpp +++ b/lib/scomplex/path_snapper.cpp @@ -1,78 +1,119 @@ -#pragma once +#include -#include +#include #include +#include +#include +#include namespace gsimp { -#include -#include -#include - struct path_snapper::impl { - tree_t point_tree; // defined in nn_utils.hpp - Graph vertex_graph; // defined in graph_utils.hpp - simplicial::SimplicialComplex* s_comp; - bool owner; + tree_t point_tree; // defined in nn_utils.hpp + graph_t vertex_graph; // defined in graph_utils.hpp + std::shared_ptr s_comp; + + impl(simplicial_complex& sc) { + s_comp.reset(&sc); + auto points = s_comp->get_points(); + make_tree(point_tree, points); + vertex_graph = calculate_one_skelleton_graph(*s_comp); + } + + impl(std::vector& pts, std::vector& cells) { + s_comp = std::shared_ptr( + new simplicial_complex(pts, cells)); + make_tree(point_tree, pts); + vertex_graph = calculate_one_skelleton_graph(*s_comp); + } - typedef std::vector, int>> geo_chain_t; - geo_chain_t get_chain_rep(std::vector path) { + ~impl(){}; + + std::vector snap_path(std::vector path) { + auto way_points = snap_points_to_indexes(point_tree, path); + auto snapped_path = complete_path(vertex_graph, way_points); + return snapped_path; + } + + std::vector> index_pairs(std::vector path) { auto vertex_path = snap_path(path); - geo_chain_t chain_rep; - for ( // - auto it = vertex_path.begin(); // - std::next(it) != vertex_path.end(); // - ++it) { - size_t s(*it), t(*(std::next(it))); - if (s < t) { - chain_rep.push_back( // - std::make_pair, int>( - std::vector{s, t}, 1)); - } else { - chain_rep.push_back( // - std::make_pair, int>( - std::vector{s, t}, -1)); - } + std::vector> pair_seq; + for (auto it = vertex_path.begin(); std::next(it) != vertex_path.end(); + ++it) { + // + // + size_t src(*it), trg(*(std::next(it))); // + if (src < trg) { // + pair_seq.push_back( // + std::make_pair, int>( // + std::vector{src, trg}, 1)); // + } else { // + pair_seq.push_back( // + std::make_pair, int>( // + std::vector{trg, src}, -1)); // + } // } - return chain_rep; + return pair_seq; } +}; + +path_snapper::path_snapper(simplicial_complex& sc) { + p_impl = std::shared_ptr(new impl(sc)); } -public : path_snapper(simplicial::SimplicialComplex& sc) { - owner = false; - s_comp = ≻ - make_tree(point_tree, s_comp->points); - vertex_graph = calculate_one_skelleton_graph(*s_comp); +path_snapper::path_snapper(std::vector& pts, + std::vector& cells) { + p_impl = std::shared_ptr(new impl(pts, cells)); } -path_snapper(std::vector& pts, std::vector& cells) { - owner = true; - s_comp = new simplicial::SimplicialComplex(pts, cells); - make_tree(point_tree, s_comp->points); - vertex_graph = calculate_one_skelleton_graph(*s_comp); +path_snapper::~path_snapper() {} + +path_snapper::path_snapper(path_snapper& other) { p_impl = other.p_impl; } + +path_snapper& path_snapper::operator=(const path_snapper& other) { + p_impl = other.p_impl; + return *this; } -~path_snapper() { - if (owner) delete s_comp; +std::vector path_snapper::snap_path_to_indices( + std::vector path) { + return p_impl->snap_path(path); } -std::vector snap_path(std::vector path) { - std::vector way_points = snap_points_to_indexes(point_tree, path); - std::vector snapped_path = complete_path(vertex_graph, way_points); - return snapped_path; +std::vector path_snapper::snap_path_to_points( + std::vector path) { + auto index_path = p_impl->snap_path(path); + std::vector point_path; + for (size_t p : index_path) + point_path.push_back(p_impl->s_comp->get_points().at(p)); + return point_path; } -vector_t get_chain_vector(std::vector path) { - auto chain_rep = get_chain_rep(path); - vector_t vector_rep(s_comp->get_dimension(1)); - for (auto chain_el : chain_rep) { - simplex_t simp = std::get<0>(chain_el); - int val = std::get<1>(chain_el); - size_t index = s_comp->get_simplex_index(simp); - vector_rep.coeffRef(index) = val; +chain_t path_snapper::get_chain(std::vector path) { + auto pair_seq = p_impl->index_pairs(path); + chain_t rep = p_impl->s_comp->new_chain(1); + for (auto pair : pair_seq) { + size_t index = p_impl->s_comp->cell_to_index(std::get<0>(pair)); + chain_val(rep, index) = std::get<1>(pair); } - return vector_rep; + return rep; } + +vector_t get_chain_vector(std::vector); +}; +/* + chain_t snap_path_to_chain(std::vector); + // interconversion + std::vector index_sequence_to_point(std::vector); + std::vector point_sequence_to_index(std::vector); + chain_t index_sequence_to_chain(std::vector); + chain_t point_sequence_to_chain(std::vector); +*/ +/* + + }; } ; + +*/ diff --git a/lib/scomplex/path_snapper.hpp b/lib/scomplex/path_snapper.hpp index 0c8cee0..90719fb 100644 --- a/lib/scomplex/path_snapper.hpp +++ b/lib/scomplex/path_snapper.hpp @@ -1,14 +1,10 @@ #pragma once -#include +#include #include namespace gsimp { -#include -#include -#include - class path_snapper { private: struct impl; @@ -17,14 +13,14 @@ class path_snapper { public: // constructors and such path_snapper(std::vector&, std::vector&); - path_snapper(simplicial::simplicial_complex&); + path_snapper(simplicial_complex&); ~path_snapper(); path_snapper(path_snapper&); - path_snapper& operator=(path_snapper&); + path_snapper& operator=(const path_snapper&); // the things we want to do std::vector snap_path_to_points(std::vector); std::vector snap_path_to_indices(std::vector); - vector_t get_chain_vector(std::vector); + chain_t get_chain(std::vector); chain_t snap_path_to_chain(std::vector); // interconversion std::vector index_sequence_to_point(std::vector); diff --git a/lib/scomplex/simplicial_complex.cpp b/lib/scomplex/simplicial_complex.cpp index f57b4b7..14ce7c1 100644 --- a/lib/scomplex/simplicial_complex.cpp +++ b/lib/scomplex/simplicial_complex.cpp @@ -1,6 +1,5 @@ #include #include -#include #include // for debuging purposes @@ -86,6 +85,18 @@ struct simplicial_complex::impl { return pow(-1, orient); } + std::vector dedupe_vec(std::vector& vec) { + std::set no_reps; + std::vector no_reps_list; + for (size_t el : vec) { + no_reps.insert(el); + } + for (size_t el : no_reps) { + no_reps_list.push_back(el); + } + return no_reps_list; + } + void calculate_matrices() { boundary_matrices = std::vector(); for (int k = 0; k < simplices.dimension(); k++) { @@ -145,6 +156,7 @@ std::vector> simplicial_complex::get_bdry_and_ind( return boundary_and_indices; }; + std::vector> simplicial_complex::get_bdry_and_ind_index( int d, size_t cell) { impl::simp_handle simp = p_impl->index_to_handle(d, cell); @@ -277,6 +289,12 @@ std::vector simplicial_complex::get_cofaces(cell_t face) { return s_cofaces; } +chain_t simplicial_complex::new_chain(int d) { + vector_t v(get_level_size(d)); + return chain_t(d,v); +} + + // TODO: get top dimensional cells from the simplicial complex to add to the // simplex tree (instead of all the simplices). simplicial_complex simplicial_complex::quotient(int char_fun(point_t)) { diff --git a/lib/scomplex/simplicial_complex.hpp b/lib/scomplex/simplicial_complex.hpp index 1fa21af..a292cc5 100644 --- a/lib/scomplex/simplicial_complex.hpp +++ b/lib/scomplex/simplicial_complex.hpp @@ -24,6 +24,8 @@ class simplicial_complex { int dimension(); bool is_quotient(); // level-wise info + chain_t new_chain(int d); + // create chains std::vector get_level(int); int get_level_size(int); // inclusion index diff --git a/lib/scomplex/types.hpp b/lib/scomplex/types.hpp index 511b445..3c0dfc9 100644 --- a/lib/scomplex/types.hpp +++ b/lib/scomplex/types.hpp @@ -1,7 +1,9 @@ +#pragma once #include #include #include + namespace gsimp { // geometric types typedef typename std::vector point_t; @@ -13,6 +15,11 @@ typedef typename Eigen::SparseVector vector_t; // chains typedef typename std::pair chain_t; +int& chain_dim(chain_t& p) { return std::get<0>(p); } +vector_t& chain_rep(chain_t& p) { return std::get<1>(p); } +double& chain_val(chain_t& p, size_t i) { return std::get<1>(p).coeffRef(i); } + + // characteristic functions // typedef std::function char_fun_t; From ed5f4d34ca15d3fb860de6785c5604f2b758b0be Mon Sep 17 00:00:00 2001 From: crvs Date: Fri, 17 Feb 2017 12:15:20 +0100 Subject: [PATCH 16/24] updated nn_utils to comply with the refactoring made --- lib/scomplex/nn_utils.hpp | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/lib/scomplex/nn_utils.hpp b/lib/scomplex/nn_utils.hpp index 756b4b3..9dd8676 100644 --- a/lib/scomplex/nn_utils.hpp +++ b/lib/scomplex/nn_utils.hpp @@ -1,18 +1,31 @@ #pragma once +#include + #include #include +// TODO: go through the includes and find out which ones are needed #include #include #include #include #include #include - #include -typedef std::vector point_t; +namespace gsimp { + +/* +types used outside: + tree_t +functions provided: + make_tree(tree_t& tree, std::vector& point_list); + nearest_neighbor(tree_t& tree, point_t point); + nearest_neighbor_index(tree_t& tree, point_t point); + snap_points(tree_t& tree, std::vector points); + snap_points_to_indexes( tree_t& tree, std::vector points); +*/ // kd-tree point type, d-dimensional point typedef CGAL::Cartesian_d::Point_d point_d; @@ -82,3 +95,4 @@ std::vector snap_points_to_indexes( // } return snapped_points; } +}; From 3cdf69f33e016e877796b81845aa204824b38ccb Mon Sep 17 00:00:00 2001 From: crvs Date: Fri, 17 Feb 2017 13:58:12 +0100 Subject: [PATCH 17/24] adapted chain calc to the current state --- CMakeLists.txt | 1 + .../{chain-calc.hpp => chain_calc.hpp} | 33 +++++++++---------- src/impl_test.cpp | 4 ++- 3 files changed, 20 insertions(+), 18 deletions(-) rename lib/scomplex/{chain-calc.hpp => chain_calc.hpp} (86%) diff --git a/CMakeLists.txt b/CMakeLists.txt index f9577bb..e66d072 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,6 +11,7 @@ project(simplicial CXX) # * CGAL under /usr/include/CGAL # set(EIGEN3_INCLUDE_DIR "/usr/local/include/eigen3") +set(EIGEN3_INCLUDE_DIR "/home/crvs/eigen-3.3.2") include_directories(${EIGEN3_INCLUDE_DIR}) set(CGAL_INCLUDE_DIR "/usr/local/include") include_directories(${CGAL_INCLUDE_DIR}) diff --git a/lib/scomplex/chain-calc.hpp b/lib/scomplex/chain_calc.hpp similarity index 86% rename from lib/scomplex/chain-calc.hpp rename to lib/scomplex/chain_calc.hpp index dd028c4..bc4677e 100644 --- a/lib/scomplex/chain-calc.hpp +++ b/lib/scomplex/chain_calc.hpp @@ -1,6 +1,8 @@ #pragma once -#include +#include +#include + #include #include #include @@ -13,27 +15,24 @@ #include #include -typedef Eigen::SparseVector vector_t; -typedef std::pair chain_t; -typedef Eigen::SparseMatrix matrix_t; -// typedef std::vector cell_t; +namespace gsimp { + +class NonZeroChain : public std::exception {}; class bounding_chain { std::vector boundary_matrices; - simplicial::SimplicialComplex* s_comp; + simplicial_complex* s_comp; bool owner; void populate_matrices(); public: - bounding_chain(simplicial::SimplicialComplex& sc); + bounding_chain(simplicial_complex& sc); bounding_chain(std::vector& points, std::vector& tris); ~bounding_chain(); chain_t get_bounding_chain(chain_t&); void draw_chain(std::string, chain_t); }; -class NonZeroChain : public std::exception {}; - //--------------------------------------- // implementation @@ -45,7 +44,7 @@ bounding_chain::~bounding_chain() { delete this->boundary_matrices.at(i); } -bounding_chain::bounding_chain(simplicial::SimplicialComplex& sc) { +bounding_chain::bounding_chain(simplicial_complex& sc) { this->owner = false; this->s_comp = ≻ this->populate_matrices(); @@ -54,15 +53,14 @@ bounding_chain::bounding_chain(simplicial::SimplicialComplex& sc) { bounding_chain::bounding_chain(std::vector& points, std::vector& tris) { this->owner = true; - this->s_comp = new simplicial::SimplicialComplex(points, tris); + this->s_comp = new simplicial_complex(points, tris); this->populate_matrices(); } void bounding_chain::populate_matrices() { for (int d = 0; d < this->s_comp->dimension(); ++d) { - auto my_matrix = this->s_comp->get_boundary(d); - this->boundary_matrices.push_back(new matrix_t(my_matrix)); - matrix_t* put_matrix = this->boundary_matrices.at(d); + auto level_matrix = this->s_comp->get_boundary_matrix(d); + this->boundary_matrices.push_back(new matrix_t(level_matrix)); } } // @@ -100,7 +98,7 @@ void bounding_chain::draw_chain(std::string filename, chain_t chain) { for (auto tri : this->s_comp->get_level(2)) { std::vector polygon_vec; for (auto k : tri) { - point_t pt = this->s_comp->points.at(k); + point_t pt = this->s_comp->get_points().at(k); point_type pt_g(pt[0], pt[1]); polygon_vec.push_back(pt_g); } @@ -141,7 +139,7 @@ void bounding_chain::draw_chain(std::string filename, chain_t chain) { if (chain_v.coeffRef(i) != 0) { std::vector line_vec; for (auto k : lines.at(i)) { - point_t pt = this->s_comp->points.at(k); + point_t pt = this->s_comp->get_points().at(k); point_type pt_g(pt[0], pt[1]); line_vec.push_back(pt_g); } @@ -180,6 +178,7 @@ chain_t bounding_chain::get_bounding_chain(chain_t& chain) { std::cout << '\n'; std::cout << result << '\n'; std::cout << chain_v << '\n'; - // throw NonZeroChain(); + throw NonZeroChain(); } } +}; diff --git a/src/impl_test.cpp b/src/impl_test.cpp index 6e82b58..4804d49 100644 --- a/src/impl_test.cpp +++ b/src/impl_test.cpp @@ -1,10 +1,12 @@ #include #include -#include #include #include #include +#include +#include +#include int my_char(gsimp::point_t point) { return 1; } From 932f6644c184374ea3cb863858e75f93fd75837b Mon Sep 17 00:00:00 2001 From: crvs Date: Fri, 17 Feb 2017 17:53:11 +0100 Subject: [PATCH 18/24] solved something with types --- lib/scomplex/types.hpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/lib/scomplex/types.hpp b/lib/scomplex/types.hpp index 3c0dfc9..590267a 100644 --- a/lib/scomplex/types.hpp +++ b/lib/scomplex/types.hpp @@ -1,4 +1,5 @@ #pragma once + #include #include #include @@ -18,9 +19,4 @@ typedef typename std::pair chain_t; int& chain_dim(chain_t& p) { return std::get<0>(p); } vector_t& chain_rep(chain_t& p) { return std::get<1>(p); } double& chain_val(chain_t& p, size_t i) { return std::get<1>(p).coeffRef(i); } - - - -// characteristic functions -// typedef std::function char_fun_t; -}; +}; // namespace gsimp From 9e7b79776f0f73f1006e56302d7fa89353bd4b15 Mon Sep 17 00:00:00 2001 From: crvs Date: Fri, 17 Feb 2017 18:37:56 +0100 Subject: [PATCH 19/24] cleaned up basic test --- src/test0.cpp | 60 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 src/test0.cpp diff --git a/src/test0.cpp b/src/test0.cpp new file mode 100644 index 0000000..5d9aea3 --- /dev/null +++ b/src/test0.cpp @@ -0,0 +1,60 @@ +/** + * SIMPLICIAL COMPLEX TEST: + * takes as argument a path to a file containing the description of a + * triangulation as output by qhull (the last column of the point coordinates + * is ignored). + * constantly asks the user to input a number which will be interpreted as a + * 1-cell index, and outputs the cofaces of the given cell, (together with the + * boundary index) and the faces of the cell. + **/ +#include +#include + +#include +#include +#include + +int my_char(gsimp::point_t point) { return 1; } + +int& do_this(int& a) { return a; } + +int main(int argc, char* argv[]) { + std::vector points_v; + std::vector cells_v; + + std::string filename; + if (argc == 1) { + std::cout << "insert file name: "; + std::cin >> filename; + } else + filename = argv[1]; + std::tie(points_v, cells_v) = parse_qhull_file(filename); + + gsimp::simplicial_complex s_comp(points_v, cells_v); + + gsimp::simplicial_complex q_comp = s_comp.quotient(my_char); + + std::cout << "complex dimension: " << s_comp.dimension() << "\n"; + for (int i = 0; i <= s_comp.dimension(); ++i) { + std::cout << "level " << i << " has " << s_comp.get_level_size(i) + << " cells\n"; + } + + size_t i; + std::cout << "get 1-cell number: "; + while (std::cin >> i) { + std::vector> my_cofaces = + s_comp.get_cof_and_ind_index(1, i); + for (auto s : my_cofaces) { + std::cout << std::get<0>(s) << ", "; + std::cout << std::get<1>(s) << " ; "; + } + std::cout << '\n'; + std::vector faces = s_comp.cell_boundary_index(1, i); + std::copy(faces.begin(), faces.end(), + std::ostream_iterator(std::cout, " ")); + std::cout << "\nget 1-cell number: "; + } + + return 0; +} From 20ebaad6be48f609dabfb83753f8235e7f9c74cb Mon Sep 17 00:00:00 2001 From: crvs Date: Fri, 17 Feb 2017 20:44:39 +0100 Subject: [PATCH 20/24] correcting errors in path snapper - added a test to the path snapping - corrected errors in the shortest path computation - added path_snapper constructor using simplicial_complex shared_poniters to avoid duplicating too much data. --- CMakeLists.txt | 7 ++-- lib/scomplex/graph_utils.hpp | 1 + lib/scomplex/path_snapper.cpp | 14 +++++++- lib/scomplex/path_snapper.hpp | 1 + src/graphtest.cpp | 15 -------- src/impl_test.cpp | 54 ---------------------------- src/libtest.cpp | 62 -------------------------------- src/main.cpp | 68 ----------------------------------- src/nn_test.cpp | 4 --- src/test0.cpp | 2 -- src/test1.cpp | 40 +++++++++++++++++++++ 11 files changed, 60 insertions(+), 208 deletions(-) delete mode 100644 src/graphtest.cpp delete mode 100644 src/impl_test.cpp delete mode 100644 src/libtest.cpp delete mode 100644 src/main.cpp delete mode 100644 src/nn_test.cpp create mode 100644 src/test1.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index e66d072..be0469b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -53,8 +53,11 @@ add_library(scomplex SHARED "./lib/scomplex/simplicial_complex.cpp") add_library(coeffflow SHARED "./lib/scomplex/coeff_flow.cpp") add_library(pathsnap SHARED "./lib/scomplex/path_snapper.cpp") -add_executable(impl_test "src/impl_test.cpp") -target_link_libraries(impl_test scomplex) +add_executable(test0 "src/test0.cpp") +target_link_libraries(test0 scomplex) + +add_executable(test1 "src/test1.cpp") +target_link_libraries(test1 scomplex pathsnap) # set(STUPID_SOURCES "./src/main.cpp") # add_executable(my_stupid_test ${STUPID_SOURCES}) diff --git a/lib/scomplex/graph_utils.hpp b/lib/scomplex/graph_utils.hpp index becf276..5ce5c26 100644 --- a/lib/scomplex/graph_utils.hpp +++ b/lib/scomplex/graph_utils.hpp @@ -81,6 +81,7 @@ std::vector shortest_path(const graph_t& g, vertex_t s, vertex_t t) { std::vector s_t_path{}; do { s_t_path.push_back(it); + it = predecessors[it]; } while (it != s); std::reverse(s_t_path.begin(), s_t_path.end()); return s_t_path; diff --git a/lib/scomplex/path_snapper.cpp b/lib/scomplex/path_snapper.cpp index 1a62637..eeb0e37 100644 --- a/lib/scomplex/path_snapper.cpp +++ b/lib/scomplex/path_snapper.cpp @@ -1,4 +1,5 @@ #include +#include #include #include @@ -13,8 +14,15 @@ struct path_snapper::impl { graph_t vertex_graph; // defined in graph_utils.hpp std::shared_ptr s_comp; + impl(std::shared_ptr p_sc) { + s_comp = p_sc; + auto points = s_comp->get_points(); + make_tree(point_tree, points); + vertex_graph = calculate_one_skelleton_graph(*s_comp); + }; + impl(simplicial_complex& sc) { - s_comp.reset(&sc); + s_comp.reset(new simplicial_complex(sc)); auto points = s_comp->get_points(); make_tree(point_tree, points); vertex_graph = calculate_one_skelleton_graph(*s_comp); @@ -57,6 +65,10 @@ struct path_snapper::impl { } }; +path_snapper::path_snapper(std::shared_ptr sc) { + p_impl = std::shared_ptr(new impl(sc)); +}; + path_snapper::path_snapper(simplicial_complex& sc) { p_impl = std::shared_ptr(new impl(sc)); } diff --git a/lib/scomplex/path_snapper.hpp b/lib/scomplex/path_snapper.hpp index 90719fb..91673ea 100644 --- a/lib/scomplex/path_snapper.hpp +++ b/lib/scomplex/path_snapper.hpp @@ -13,6 +13,7 @@ class path_snapper { public: // constructors and such path_snapper(std::vector&, std::vector&); + path_snapper(std::shared_ptr); path_snapper(simplicial_complex&); ~path_snapper(); path_snapper(path_snapper&); diff --git a/src/graphtest.cpp b/src/graphtest.cpp deleted file mode 100644 index 6329c70..0000000 --- a/src/graphtest.cpp +++ /dev/null @@ -1,15 +0,0 @@ - -#include - -#include "boost/graph/adjacency_list.hpp" -#include "boost/graph/dijkstra_shortest_paths.hpp" - -typedef boost::adjacency_list - Graph; - -int main() { - Graph g(10); - boost::add_edge(0, 1, g); - - return 0; -} diff --git a/src/impl_test.cpp b/src/impl_test.cpp deleted file mode 100644 index 4804d49..0000000 --- a/src/impl_test.cpp +++ /dev/null @@ -1,54 +0,0 @@ -#include -#include - -#include -#include -#include -#include -#include -#include - -int my_char(gsimp::point_t point) { return 1; } - -int& do_this(int& a) { return a; } - -int main(int argc, char* argv[]) { - std::vector points_v; - std::vector cells_v; - - std::string filename; - if (argc == 1) { - std::cout << "insert file name: "; - std::cin >> filename; - } else - filename = argv[1]; - std::tie(points_v, cells_v) = parse_qhull_file(filename); - - gsimp::simplicial_complex s_comp(points_v, cells_v); - - gsimp::simplicial_complex q_comp = s_comp.quotient(my_char); - - std::cout << "complex dimension: " << s_comp.dimension() << "\n"; - for (int i = 0; i <= s_comp.dimension(); ++i) { - std::cout << "level " << i << " has " << s_comp.get_level_size(i) - << " cells\n"; - } - - size_t i; - std::cout << "get 1-cell number: "; - while (std::cin >> i) { - std::vector> my_cofaces = - s_comp.get_cof_and_ind_index(2, i); - for (auto s : my_cofaces) { - std::cout << std::get<0>(s) << ", "; - std::cout << std::get<1>(s) << " ,, "; - } - std::cout << '\n'; - std::vector faces = s_comp.cell_boundary_index(2, i); - std::copy(faces.begin(), faces.end(), - std::ostream_iterator(std::cout, " ")); - std::cout << "\nget 1-cell number: "; - } - - return 0; -} diff --git a/src/libtest.cpp b/src/libtest.cpp deleted file mode 100644 index c9ca754..0000000 --- a/src/libtest.cpp +++ /dev/null @@ -1,62 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include "boost/graph/graph_traits.hpp" - -#include -#include -#include -#include - -typedef std::vector point_t; -typedef std::vector cell_t; - -int f(point_t a) { return a.at(0) < .5 ? 0 : 1; } - -int main() { - std::vector points_v; - std::vector cells_v; - - std::string filename; - std::cout << "insert file name: "; - std::cin >> filename; - std::tie(points_v, cells_v) = parse_qhull_file(filename); - - simplicial::SimplicialComplex rsc(points_v, cells_v); - auto qsc(rsc.quotient(f)); - snap::path_snapper snapper(rsc); - - auto path = snapper.snap_path({{-1, -1}, {1, 0}, {1, 1}}); - std::copy(path.begin(), path.end(), - std::ostream_iterator(std::cout, " ")); - std::cout << '\n'; - - bounding_chain chain_calc(points_v, cells_v); - - path = snapper.snap_path({{-1, -1}, {1, 0}, {1, 1}, {-1, -1}}); - std::copy(path.begin(), path.end(), - std::ostream_iterator(std::cout, " ")); - std::cout << '\n'; - vector_t path_v( - snapper.get_chain_vector({{-1, -1}, {1, 0}, {1, 1}, {-1, -1}})); - - // try { - chain_t chain(2, path_v); - chain = chain_calc.get_bounding_chain(chain); - std::cout << path_v << '\n'; - vector_t b_chain = round_vec(std::get<1>(chain)); - std::cout << b_chain; - // std::cout << chain_r << '\n'; - std::cout << round_vec(rsc.get_boundary(1) * std::get<1>(chain)) << '\n'; - std::cout << round_vec(rsc.get_boundary(1) * b_chain) << '\n'; - - //} catch (NonZeroChain) { - std::cout << "not a zero chain"; - //} - - return 0; -} diff --git a/src/main.cpp b/src/main.cpp deleted file mode 100644 index c296cbb..0000000 --- a/src/main.cpp +++ /dev/null @@ -1,68 +0,0 @@ -#include - -#include -#include - -int main() { - struct Myoptions : Gudhi::Simplex_tree_options_full_featured { - typedef int Vertex_handle; - }; - Gudhi::Simplex_tree simplexTree; - auto d = {0, 1, 2}; - simplexTree.insert_simplex_and_subfaces(d); - - auto range = simplexTree.complex_simplex_range(); - int count = 0; - - for (auto s : range) { - simplexTree.assign_key(s, count++); - auto v_range = simplexTree.simplex_vertex_range(s); - std::cout << simplexTree.filtration(s) << " , " << simplexTree.key(s) - << " : "; - for (auto v : v_range) { - std::cout << v << " "; - } - std::cout << ": "; - auto sbd = simplexTree.boundary_simplex_range(s); - for (auto sb : sbd) { - auto v_range = simplexTree.simplex_vertex_range(sb); - for (auto v : v_range) { - std::cout << v << " "; - } - std::cout << " ; "; - } - std::cout << std::endl; - } - - // testing the effects of reinsertion - simplexTree.insert_simplex_and_subfaces(d); - for (auto s : range) { - auto v_range = simplexTree.simplex_vertex_range(s); - std::cout << simplexTree.filtration(s) << " , " << simplexTree.key(s) - << " : "; - for (auto v : v_range) { - std::cout << v << " "; - } - std::cout << ": "; - auto sbd = simplexTree.boundary_simplex_range(s); - for (auto sb : sbd) { - auto v_range = simplexTree.simplex_vertex_range(sb); - for (auto v : v_range) { - std::cout << v << " "; - } - std::cout << " ; "; - } - std::cout << std::endl; - } - - auto p = simplexTree.find({0, 1, 2}); - for (auto v : simplexTree.simplex_vertex_range(p)) { - std::cout << v << " "; - } - std::cout << std::endl; - - std::vector vec; - std::cout << "vector size " << vec.size() << std::endl; - - return 0; -} diff --git a/src/nn_test.cpp b/src/nn_test.cpp deleted file mode 100644 index f6fb1f7..0000000 --- a/src/nn_test.cpp +++ /dev/null @@ -1,4 +0,0 @@ -#include -#include - -int main() { return 0; } diff --git a/src/test0.cpp b/src/test0.cpp index 5d9aea3..8637fd8 100644 --- a/src/test0.cpp +++ b/src/test0.cpp @@ -32,8 +32,6 @@ int main(int argc, char* argv[]) { gsimp::simplicial_complex s_comp(points_v, cells_v); - gsimp::simplicial_complex q_comp = s_comp.quotient(my_char); - std::cout << "complex dimension: " << s_comp.dimension() << "\n"; for (int i = 0; i <= s_comp.dimension(); ++i) { std::cout << "level " << i << " has " << s_comp.get_level_size(i) diff --git a/src/test1.cpp b/src/test1.cpp new file mode 100644 index 0000000..ac88e1f --- /dev/null +++ b/src/test1.cpp @@ -0,0 +1,40 @@ +/** + **/ +#include +#include + +#include +#include +#include +#include + +int my_char(gsimp::point_t point) { return 1; } + +int& do_this(int& a) { return a; } + +int main(int argc, char* argv[]) { + std::vector points_v; + std::vector cells_v; + + std::string filename; + if (argc == 1) { + std::cout << "insert file name: "; + std::cin >> filename; + } else + filename = argv[1]; + std::tie(points_v, cells_v) = parse_qhull_file(filename); + + auto s_comp=std::shared_ptr(new gsimp::simplicial_complex(points_v, cells_v)); + + std::cout << "complex dimension: " << s_comp->dimension() << "\n"; + for (int i = 0; i <= s_comp->dimension(); ++i) { + std::cout << "level " << i << " has " << s_comp->get_level_size(i) + << " cells\n"; + } + + gsimp::path_snapper p_snap(s_comp); + std::vector snapped = p_snap.snap_path_to_indices({{-2,-2},{-2,2},{2,2},{-2,-2}}); + std::copy(snapped.begin(),snapped.end(),std::ostream_iterator(std::cout," ")); + + return 0; +} From 86807e42bb3a6a3e5012849c26b15ec725bd0f17 Mon Sep 17 00:00:00 2001 From: crvs Date: Fri, 17 Feb 2017 21:13:55 +0100 Subject: [PATCH 21/24] finished implementations in path_snapper header --- lib/scomplex/path_snapper.cpp | 49 +++++++++++++++++++++-------------- lib/scomplex/path_snapper.hpp | 1 - src/test1.cpp | 12 ++++++--- 3 files changed, 38 insertions(+), 24 deletions(-) diff --git a/lib/scomplex/path_snapper.cpp b/lib/scomplex/path_snapper.cpp index eeb0e37..23ce003 100644 --- a/lib/scomplex/path_snapper.cpp +++ b/lib/scomplex/path_snapper.cpp @@ -87,11 +87,6 @@ path_snapper& path_snapper::operator=(const path_snapper& other) { return *this; } -std::vector path_snapper::snap_path_to_indices( - std::vector path) { - return p_impl->snap_path(path); -} - std::vector path_snapper::snap_path_to_points( std::vector path) { auto index_path = p_impl->snap_path(path); @@ -101,7 +96,12 @@ std::vector path_snapper::snap_path_to_points( return point_path; } -chain_t path_snapper::get_chain(std::vector path) { +std::vector path_snapper::snap_path_to_indices( + std::vector path) { + return p_impl->snap_path(path); +} + +chain_t path_snapper::snap_path_to_chain(std::vector path) { auto pair_seq = p_impl->index_pairs(path); chain_t rep = p_impl->s_comp->new_chain(1); for (auto pair : pair_seq) { @@ -111,21 +111,30 @@ chain_t path_snapper::get_chain(std::vector path) { return rep; } -vector_t get_chain_vector(std::vector); -}; -/* - chain_t snap_path_to_chain(std::vector); - // interconversion - std::vector index_sequence_to_point(std::vector); - std::vector point_sequence_to_index(std::vector); - chain_t index_sequence_to_chain(std::vector); - chain_t point_sequence_to_chain(std::vector); -*/ -/* +std::vector path_snapper::index_sequence_to_point( + std::vector ind_path) { + std::vector pt_path; + for (size_t ind : ind_path) + pt_path.push_back(p_impl->s_comp->get_points().at(ind)); + return pt_path; +} +std::vector path_snapper::point_sequence_to_index( + std::vector pt_path) { + std::vector ind_path = + snap_points_to_indexes(p_impl->point_tree, pt_path); + return ind_path; +} -}; +chain_t path_snapper::index_sequence_to_chain(std::vector ind_path) { + chain_t chain = p_impl->s_comp->new_chain(1); + auto it = ind_path.begin(); + for (; std::next(it) != ind_path.end(); ++it) + chain_val(chain, p_impl->s_comp->cell_to_index({*it, *std::next(it)})); + return chain; } -; -*/ +chain_t path_snapper::point_sequence_to_chain(std::vector pt_path) { + return index_sequence_to_chain(point_sequence_to_index(pt_path)); +} +}; diff --git a/lib/scomplex/path_snapper.hpp b/lib/scomplex/path_snapper.hpp index 91673ea..43025ea 100644 --- a/lib/scomplex/path_snapper.hpp +++ b/lib/scomplex/path_snapper.hpp @@ -21,7 +21,6 @@ class path_snapper { // the things we want to do std::vector snap_path_to_points(std::vector); std::vector snap_path_to_indices(std::vector); - chain_t get_chain(std::vector); chain_t snap_path_to_chain(std::vector); // interconversion std::vector index_sequence_to_point(std::vector); diff --git a/src/test1.cpp b/src/test1.cpp index ac88e1f..50c2431 100644 --- a/src/test1.cpp +++ b/src/test1.cpp @@ -24,7 +24,8 @@ int main(int argc, char* argv[]) { filename = argv[1]; std::tie(points_v, cells_v) = parse_qhull_file(filename); - auto s_comp=std::shared_ptr(new gsimp::simplicial_complex(points_v, cells_v)); + auto s_comp = std::shared_ptr( + new gsimp::simplicial_complex(points_v, cells_v)); std::cout << "complex dimension: " << s_comp->dimension() << "\n"; for (int i = 0; i <= s_comp->dimension(); ++i) { @@ -33,8 +34,13 @@ int main(int argc, char* argv[]) { } gsimp::path_snapper p_snap(s_comp); - std::vector snapped = p_snap.snap_path_to_indices({{-2,-2},{-2,2},{2,2},{-2,-2}}); - std::copy(snapped.begin(),snapped.end(),std::ostream_iterator(std::cout," ")); + std::vector snapped = + p_snap.snap_path_to_indices({{-2, -2}, {-2, 2}, {2, 2}, {-2, -2}}); + std::copy(snapped.begin(), snapped.end(), + std::ostream_iterator(std::cout, " ")); + std::cout << '\n'; + + auto chain = p_snap.index_sequence_to_chain(snapped); return 0; } From 097a2e3c7e8ed9d0ba1a4ad23f12b2a606ee90d7 Mon Sep 17 00:00:00 2001 From: crvs Date: Fri, 17 Feb 2017 21:45:46 +0100 Subject: [PATCH 22/24] removed old stuff --- lib/scomplex/Simplicial_Complex.h | 235 ------------------------------ lib/scomplex/trace.h | 55 ------- lib/scomplex/utils.h | 5 - 3 files changed, 295 deletions(-) delete mode 100644 lib/scomplex/Simplicial_Complex.h delete mode 100644 lib/scomplex/trace.h delete mode 100644 lib/scomplex/utils.h diff --git a/lib/scomplex/Simplicial_Complex.h b/lib/scomplex/Simplicial_Complex.h deleted file mode 100644 index 17a683b..0000000 --- a/lib/scomplex/Simplicial_Complex.h +++ /dev/null @@ -1,235 +0,0 @@ -#pragma once - -#include - -#include -#include - -#include - -#include -#include -#include - -typedef std::vector point_t; -namespace simplicial { - -class NoChain {}; - -class SimplicialComplex { - private: - // the type of matrix to be used - typedef typename Eigen::SparseMatrix matrix_t; - typedef typename Eigen::SparseVector vector_t; - - // simplex tree options - struct SimpleOptions : Gudhi::Simplex_tree_options_full_featured { - typedef size_t Vertex_handle; - }; - typedef typename Gudhi::Simplex_tree SimplexTree; - typedef typename Gudhi::Simplex_tree::Simplex_handle - Simplex_handle; - - typedef std::vector level_t; - typedef std::vector levels_t; - typedef std::vector simplex_t; - - public: - std::vector get_level(int dimen) { - std::vector level; - for (auto simp : simplices.complex_simplex_range()) { - if (simplices.dimension(simp) == dimen) { - simplex_t v_simp; - for (auto v : simplices.simplex_vertex_range(simp)) { - v_simp.push_back(v); - } - level.push_back(v_simp); - } - } - return level; - } - - std::vector points; - SimplexTree simplices; - - int get_dimension(int level) { return dimensions.at(level); } - - SimplicialComplex() : geometric_q(true) {} - - /** - * @brief builds a simplicial complex object out of points and "triangles" - * - * @param points_a points - * @param tris the top dimensional simplicial complex structure - */ - SimplicialComplex(std::vector points_a, - std::vector> tris) - : geometric_q(true), points(points_a) { - for (auto s : tris) { - // cant have simplices with repeated vertices so we dedupe the lists - simplices.insert_simplex_and_subfaces(dedupe_list(s)); - int d = s.size() - 1; - if (simplices.dimension() < d) { - simplices.set_dimension(d); - } - } - - // we need to count i--simplices in order to get a key for them - std::vector count = std::vector(simplices.dimension() + 1, 0); - levels_t levels; - levels = levels_t(simplices.dimension() + 1, level_t()); - - auto simplex_range = simplices.complex_simplex_range(); - for (auto s : simplex_range) { - int d = simplices.dimension(s); - simplices.assign_key(s, count.at(d)++); - levels.at(d).push_back(&s); - auto v_range = simplices.simplex_vertex_range(s); - auto sbd = simplices.boundary_simplex_range(s); - } - dimensions = count; - } - - matrix_t get_boundary(int d) { - // uninstantiated boundary matrices - if (boundary_matrices.size() == 0) { - calculate_matrices(); - } - - // now they have to be instantiated, get them - if (0 <= d && d < boundary_matrices.size()) { - return boundary_matrices.at(d); - } else { - throw NoChain(); - } - } - - SimplicialComplex quotient(int f(point_t)) { - int n_points = 1; - auto corresp = std::vector(); - - // making the correspondence between points in the original complex - // and - // points in the quotient. point index 0 corresponds to the - // "virtual" - // point (in case it exists). - std::vector q_points; - bool geometric = true; - - for (auto p : points) { - if (f(p) != 0) { - corresp.push_back(n_points++); - q_points.push_back(p); - } else { - if (geometric) { - // need to realize that it is not geometric, (flip - // geometric - // and add the virtual point - auto pos = q_points.begin(); - q_points.insert(pos, point_t()); - geometric = false; - } - corresp.push_back(0); - } - } - if (geometric) { - for (int i = 0; i < corresp.size(); i++) { - corresp.at(i) -= 1; - } - } - - // producing list of simplices - std::vector> simp_list; - - for (auto s : simplices.complex_simplex_range()) { - std::vector s_q; - for (int v : simplices.simplex_vertex_range(s)) { - s_q.push_back(corresp.at(v)); - } - simp_list.push_back(dedupe_list(s_q)); - } - - auto quotient_sc = SimplicialComplex(q_points, simp_list); - quotient_sc.geometric_q = geometric; - - return quotient_sc; - } - - /** - * @brief function deciding if a simplicial complex is geometric, note - * that - * quotient simplicial complexes are not geometric. - */ - bool is_geometric() { return geometric_q; } - - size_t get_simplex_index(simplex_t simp) { - auto sh = simplices.find(simp); - return simplices.key(sh); - } - - int dimension() { return simplices.dimension(); } - - private: - bool geometric_q; - std::vector dimensions; - - std::vector boundary_matrices; - - std::set vertex_set(Simplex_handle s) { - std::set v_set; - for (auto v : simplices.simplex_vertex_range(s)) { - v_set.insert(v); - } - return v_set; - } - - // calculate the index of s_1 in the boundary of s_2 - int boundary_index(Simplex_handle s_1, Simplex_handle s_2) { - auto it_1 = simplices.simplex_vertex_range(s_1).begin(); - auto end_1 = simplices.simplex_vertex_range(s_1).end(); - // - auto it_2 = simplices.simplex_vertex_range(s_2).begin(); - auto end_2 = simplices.simplex_vertex_range(s_2).end(); - // - int orient = 0; - int unmatch = 0; - // - for (int p = 0; p <= simplices.dimension(s_2); p++) { - if (*it_1 != *it_2) { - orient = p++; - unmatch++; - } else { - if (it_1 != end_1) { - it_1++; - } - it_2++; - } - } - return pow(-1, orient); - } - - /** - * @brief calculate the boundary matrices - */ - void calculate_matrices() { - boundary_matrices = std::vector(); - for (int k = 0; k < simplices.dimension(); k++) { - boundary_matrices.push_back( - matrix_t(dimensions.at(k), dimensions.at(k + 1))); - } - for (auto s : simplices.complex_simplex_range()) { - int j = simplices.key(s); - for (auto bs : simplices.boundary_simplex_range(s)) { - int i = simplices.key(bs); - int k = simplices.dimension(bs); - boundary_matrices.at(k).coeffRef(i, j) = boundary_index(bs, s); - } - } - } - - bool is_point(point_t pt) { - return (geometric_q ? true : not(pt.size() == 0)); - } - -}; // class SimplicialComplex -}; // namespace simplicial diff --git a/lib/scomplex/trace.h b/lib/scomplex/trace.h deleted file mode 100644 index e747213..0000000 --- a/lib/scomplex/trace.h +++ /dev/null @@ -1,55 +0,0 @@ -// here we code the trace_ objects they keep a path -#ifndef TRACE_H -#define TRACE_H - -#include -#include -#include - -namespace trace { -typedef std::vector point_t; -class trace { - private: - size_t dimension_; - size_t length_; - std::vector trace_; - - public: - trace(std::vector& list) { - trace_ = std::vector(list); - - // error if we make an empty trace - assert(trace_.size() > 0); - - dimension_ = trace_.at(0).size(); - length_ = trace_.size(); - } - - trace(trace& other) { - trace_ = std::vector(other.trace_); - dimension_ = trace_.at(0).size(); - length_ = trace_.size(); - } - - size_t dimension() { return dimension_; }; - size_t size() { return length_; } - std::vector get_trace() { return trace_; } - - std::string print() { - std::string out_str; - out_str = "dim: "; - out_str += std::to_string(dimension_); - out_str += " size: "; - out_str += std::to_string(length_); - out_str += " points"; - for (point_t pt : trace_) { - out_str += ":"; - for (double coord : pt) { - out_str += std::to_string(coord) + " "; - } - } - out_str += std::to_string(length_); - } -}; -}; -#endif diff --git a/lib/scomplex/utils.h b/lib/scomplex/utils.h deleted file mode 100644 index d34a5ec..0000000 --- a/lib/scomplex/utils.h +++ /dev/null @@ -1,5 +0,0 @@ -#pragma once - -#include "scomplex/Simplicial_Complex.h" -#include "scomplex/graph_utils.hpp" -#include "scomplex/nn_utils.hpp" From ab0d3540d70f81f7a0b5f721d8bb66a61ecc116f Mon Sep 17 00:00:00 2001 From: crvs Date: Fri, 17 Feb 2017 21:47:05 +0100 Subject: [PATCH 23/24] more cleanup --- lib/scomplex/path-snapping.hpp | 28 ---------------------------- 1 file changed, 28 deletions(-) delete mode 100644 lib/scomplex/path-snapping.hpp diff --git a/lib/scomplex/path-snapping.hpp b/lib/scomplex/path-snapping.hpp deleted file mode 100644 index ca5bce1..0000000 --- a/lib/scomplex/path-snapping.hpp +++ /dev/null @@ -1,28 +0,0 @@ -#pragma once - -namespace snap { - -#include -#include -#include -#include -#include - -typedef Eigen::SparseVector vector_t; - -class path_snapper { - private: - struct impl; - std::shared_ptr p_impl; - - public: - path_snapper(std::vector& pts, std::vector& cells); - path_snapper(simplicial::simplicial_complex&); - ~path_snapper(); - path_snapper(path_snapper& other); - path_snapper& operator=(path_snapper& other); - - std::vector snap_path(std::vector path); - vector_t get_chain_vector(std::vector path); -}; -}; From 5111aedae3d32f768596d35f9bcfae29fcc79aea Mon Sep 17 00:00:00 2001 From: crvs Date: Sat, 25 Feb 2017 00:58:39 +0100 Subject: [PATCH 24/24] finished coeff_flow as well as pictures --- CMakeLists.txt | 12 ++ img-test/Makefile | 5 - img-test/test0.cpp | 195 ----------------- lib/scomplex/chain_calc.hpp | 197 +++++++++--------- .../{coeff_flow.cpp => coeff_flow.hpp} | 53 +++-- lib/scomplex/path_snapper.cpp | 10 +- lib/scomplex/simplicial_complex.cpp | 2 + lib/scomplex/types.hpp | 49 ++++- src/test0.cpp | 7 +- src/test1.cpp | 9 + src/test2.cpp | 60 ++++++ src/test3.cpp | 48 +++++ src/test4.cpp | 61 ++++++ src/test5.cpp | 58 ++++++ src/utils.hpp | 2 + 15 files changed, 439 insertions(+), 329 deletions(-) delete mode 100644 img-test/Makefile delete mode 100644 img-test/test0.cpp rename lib/scomplex/{coeff_flow.cpp => coeff_flow.hpp} (66%) create mode 100644 src/test2.cpp create mode 100644 src/test3.cpp create mode 100644 src/test4.cpp create mode 100644 src/test5.cpp create mode 100644 src/utils.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index be0469b..7b732e5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -59,6 +59,18 @@ target_link_libraries(test0 scomplex) add_executable(test1 "src/test1.cpp") target_link_libraries(test1 scomplex pathsnap) +add_executable(test2 "src/test2.cpp") +target_link_libraries(test2 scomplex pathsnap) + +add_executable(test3 "src/test3.cpp") +target_link_libraries(test3 scomplex pathsnap) + +add_executable(test4 "src/test4.cpp") +target_link_libraries(test4 scomplex pathsnap) + +add_executable(test5 "src/test5.cpp") +target_link_libraries(test5 scomplex pathsnap) + # set(STUPID_SOURCES "./src/main.cpp") # add_executable(my_stupid_test ${STUPID_SOURCES}) diff --git a/img-test/Makefile b/img-test/Makefile deleted file mode 100644 index 8339643..0000000 --- a/img-test/Makefile +++ /dev/null @@ -1,5 +0,0 @@ -image: a.out - ./a.out - -a.out: test0.cpp - g++ test0.cpp -ggdb --std=c++11 -I/usr/include/eigen3 -lpng diff --git a/img-test/test0.cpp b/img-test/test0.cpp deleted file mode 100644 index cc47da1..0000000 --- a/img-test/test0.cpp +++ /dev/null @@ -1,195 +0,0 @@ -// generic image libraries -#include -#include -// vector arythmetics -#include -using namespace boost::gil; -#include -#include -#include - -class Segment { - Eigen::Vector2d p0, p1; - Eigen::Vector2d director, normal; - - public: - Segment(Eigen::Vector2d arg_p0, Eigen::Vector2d arg_p1) { - p0 = arg_p0; - p1 = arg_p1; - - assert(length() > 0 && "Points have to be different"); - - director = p1 - p0; - director = (1 / director.norm()) * director; - normal << director(1), -director(0); - normal = (1 / normal.norm()) * normal; - } - float length() { return (p1 - p0).norm(); } - - private: - bool is_in_segment(Eigen::Vector2d point, float thickness) { - bool is_in_segment_v; - is_in_segment_v = fabs(normal.dot(point - p0)) < thickness / 2 // - && director.dot(point - p0) > -thickness / 2 // - && -director.dot(point - p1) > -thickness / 2; - return is_in_segment_v; - } - - public: - void draw(rgba8_image_t& img, // - const rgba8_pixel_t color, // - const float thickness // - ) { - int width = img.width(); - int height = img.height(); - auto v = view(img); - for (int x = 0; x < width; x++) { - for (int y = 0; y < height; y++) { - Eigen::Vector2d point; - point << x, y; - if (is_in_segment(point, thickness)) { - v(x, y) = color; - } - } - } - } -}; - -class Triangle { - Eigen::Vector2d t0, t1, t2; - - private: - Eigen::Vector2d interior_0, interior_1, interior_2; - - public: - Triangle(const Eigen::Vector2d arg_t0, // - const Eigen::Vector2d arg_t1, // - const Eigen::Vector2d arg_t2) { - t0 = arg_t0; - t1 = arg_t1; - t2 = arg_t2; - - Eigen::Vector2d director_0 = (t1 - t0); - Eigen::Vector2d director_1 = (t2 - t1); - Eigen::Vector2d director_2 = (t0 - t2); - - interior_0 << director_0(1), -1 * director_0(0); - if (interior_0.dot(director_1) < 0) { - interior_0 = -1 * interior_0; - } - assert(interior_0.dot(director_1) > 0 && - "points are not in general position"); - interior_0 = (1 / interior_0.norm()) * interior_0; - - interior_1 << director_1(1), -1 * director_1(0); - if (interior_1.dot(director_2) < 0) { - interior_1 = -1 * interior_1; - } - interior_1 = (1 / interior_1.norm()) * interior_1; - - interior_2 << director_2(1), -1 * director_2(0); - if (interior_2.dot(director_0) < 0) { - interior_2 = -1 * interior_2; - } - interior_2 = (1 / interior_2.norm()) * interior_2; - } - - private: - bool is_in_border(const Eigen::Vector2d point, const float thickness) { - bool is_in_border_v = false; - // do it all at once so that the compiler can optimize - is_in_border_v = - ((fabs(interior_0.dot(point - t0)) < thickness / 2) // - && (interior_1.dot(point - t1) > -thickness / 2) // - && (interior_2.dot(point - t2) > -thickness / 2)) || // - ((fabs(interior_1.dot(point - t1)) < thickness / 2) // - && (interior_2.dot(point - t2) > -thickness / 2) // - && (interior_0.dot(point - t0) > -thickness / 2)) || // - ((fabs(interior_2.dot(point - t2)) < thickness / 2) // - && (interior_1.dot(point - t1) > -thickness / 2) // - && (interior_0.dot(point - t0) > -thickness / 2)); - return is_in_border_v; - }; - - bool is_inside(const Eigen::Vector2d point) { - bool is_inside_v = (interior_0.dot(point - t0) > 0) // - && (interior_1.dot(point - t1) > 0) // - && (interior_2.dot(point - t2) > 0); - return is_inside_v; - } - - void draw_border(rgba8_image_t& img, // - const rgba8_pixel_t color, // - const float thickness // - ) { - // getting the view - auto v = view(img); - // bounds - int width = img.width(); - int height = img.height(); - - for (int x = 0; x < width; x++) { - for (int y = 0; y < height; y++) { - Eigen::Vector2d point; - point << x, y; - if (is_in_border(point, thickness)) { - v(x, y) = color; - } - } - } - }; - - void draw_interior(rgba8_image_t& img, // - const rgba8_pixel_t color // - ) { - // getting the view - auto v = view(img); - // bounds - int width = img.width(); - int height = img.height(); - - for (int x = 0; x < width; x++) { - for (int y = 0; y < height; y++) { - Eigen::Vector2d point; - point << x, y; - if (is_inside(point)) { - v(x, y) = color; - } - } - } - } - - public: - void draw(rgba8_image_t& img, // - const rgba8_pixel_t interior_color, // - const rgba8_pixel_t border_color, // - const float border_thickness) // - { - // first draw the interior of the triangle - draw_interior(img, interior_color); - // after the interior add the border - draw_border(img, border_color, border_thickness); - }; - - float area() { - Eigen::Vector2d d1{t1 - t0}, d2{t2 - t0}; - return (d1(0) * d2(1) - d1(1) * d2(0)) / 2; - } -}; - -int main() { - Eigen::Vector2d a0(10, 10), a1(100, 10), a2(10, 100); - Triangle T(a0, a1, a2); - rgba8_image_t img(512, 512); - rgba8_pixel_t red(255, 0, 0, 128); - fill_pixels(view(img), red); - T.draw(img, {0, 255, 255, 128}, {0, 0, 0, 255}, 2); - Triangle T2({100, 10}, {10, 100}, {100, 100}); - T2.draw(img, {0, 255, 255, 128}, {0, 0, 0, 255}, 2); - - Segment s({100, 230}, {100, 230}); - std::cout << std::endl - << T.area() << std::endl; - - png_write_view("redsquare.png", const_view(img)); -} diff --git a/lib/scomplex/chain_calc.hpp b/lib/scomplex/chain_calc.hpp index bc4677e..b88f9fc 100644 --- a/lib/scomplex/chain_calc.hpp +++ b/lib/scomplex/chain_calc.hpp @@ -10,6 +10,7 @@ #include #include +#include #include #include #include @@ -17,50 +18,49 @@ namespace gsimp { -class NonZeroChain : public std::exception {}; +class non_zero_chain : public std::exception {}; +class cant_draw_dimension : public std::exception {}; class bounding_chain { - std::vector boundary_matrices; - simplicial_complex* s_comp; - bool owner; + std::shared_ptr s_comp; + std::vector> boundary_matrices; void populate_matrices(); public: bounding_chain(simplicial_complex& sc); + bounding_chain(std::shared_ptr sc); bounding_chain(std::vector& points, std::vector& tris); ~bounding_chain(); chain_t get_bounding_chain(chain_t&); - void draw_chain(std::string, chain_t); + // void draw_chain(std::string, chain_t); }; //--------------------------------------- // implementation -bounding_chain::~bounding_chain() { - if (this->owner) delete this->s_comp; +bounding_chain::~bounding_chain() {} - // we always own the matrices - for (int i = this->boundary_matrices.size() - 1; i >= 0; --i) - delete this->boundary_matrices.at(i); +bounding_chain::bounding_chain(std::shared_ptr sc) { + s_comp = sc; + populate_matrices(); } bounding_chain::bounding_chain(simplicial_complex& sc) { - this->owner = false; - this->s_comp = ≻ - this->populate_matrices(); + s_comp.reset(new simplicial_complex(sc)); + populate_matrices(); } bounding_chain::bounding_chain(std::vector& points, std::vector& tris) { - this->owner = true; - this->s_comp = new simplicial_complex(points, tris); - this->populate_matrices(); + s_comp.reset(new simplicial_complex(points, tris)); + populate_matrices(); } void bounding_chain::populate_matrices() { - for (int d = 0; d < this->s_comp->dimension(); ++d) { - auto level_matrix = this->s_comp->get_boundary_matrix(d); - this->boundary_matrices.push_back(new matrix_t(level_matrix)); + for (int d = 0; d < s_comp->dimension(); ++d) { + auto level_matrix = s_comp->get_boundary_matrix(d); + boundary_matrices.push_back( + std::unique_ptr(new matrix_t(level_matrix))); } } // @@ -85,100 +85,99 @@ bool equals(vector_t vec1, vector_t vec2) { return true; } -// TODO: drowing facilities -// - draw chain -// - draw 1-chain with arrows -// - put it into its own header file -void bounding_chain::draw_chain(std::string filename, chain_t chain) { - typedef boost::geometry::model::d2::point_xy point_type; - typedef boost::geometry::model::ring polygon_type; - - std::vector polys; - - for (auto tri : this->s_comp->get_level(2)) { - std::vector polygon_vec; - for (auto k : tri) { - point_t pt = this->s_comp->get_points().at(k); - point_type pt_g(pt[0], pt[1]); - polygon_vec.push_back(pt_g); - } - polygon_type poly(polygon_vec.begin(), polygon_vec.end()); - polys.push_back(poly); - } - - // now we want to draw this - std::ofstream svg(filename.c_str()); - - boost::geometry::svg_mapper mapper(svg, 150, 150); - - int chain_d = std::get<0>(chain); - vector_t chain_v = std::get<1>(chain); - for (int i = 0; i < chain_v.rows(); ++i) { - polygon_type poly = polys.at(i); - mapper.add(poly); - if (chain_d == 2) { - if (chain_v.coeffRef(i) == 0) - mapper.map( - poly, - "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153," - "204,0);stroke-width:2"); - else - mapper.map( - poly, - "fill-opacity:0.5;fill:rgb(230,104,58);stroke:rgb(230," - "104,58);stroke-width:2"); - } else { - mapper.map(poly, - "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153," - "204,0);stroke-width:2"); - } - - if (chain_d == 1) { - auto lines = this->s_comp->get_level(1); - for (int i = 0; i < lines.size(); ++i) { - if (chain_v.coeffRef(i) != 0) { - std::vector line_vec; - for (auto k : lines.at(i)) { - point_t pt = this->s_comp->get_points().at(k); - point_type pt_g(pt[0], pt[1]); - line_vec.push_back(pt_g); - } - polygon_type line_g(line_vec.begin(), line_vec.end()); - mapper.map( - line_g, - "fill-opacity:0.5;fill:rgb(250,0,0);stroke:rgb(250," - "0,0);stroke-width:2"); - } - } - } - } -} +// // TODO: drowing facilities +// // - draw chain +// // - draw 1-chain with arrows +// // - put it into its own header file +// void bounding_chain::draw_chain(std::string filename, chain_t chain) { +// typedef boost::geometry::model::d2::point_xy point_type; +// typedef boost::geometry::model::ring polygon_type; +// +// std::vector polys; +// +// for (auto tri : s_comp->get_level(2)) { +// std::vector polygon_vec; +// for (auto k : tri) { +// point_t pt = s_comp->get_points().at(k); +// point_type pt_g(pt[0], pt[1]); +// polygon_vec.push_back(pt_g); +// } +// polygon_type poly(polygon_vec.begin(), polygon_vec.end()); +// polys.push_back(poly); +// } +// +// // now we want to draw this +// std::ofstream svg(filename.c_str()); +// +// boost::geometry::svg_mapper mapper(svg, 150, 150); +// +// int chain_d = std::get<0>(chain); +// vector_t chain_v = std::get<1>(chain); +// for (int i = 0; i < chain_v.rows(); ++i) { +// polygon_type poly = polys.at(i); +// mapper.add(poly); +// if (chain_d == 2) { +// if (chain_v.coeffRef(i) == 0) +// mapper.map( +// poly, +// "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153," +// "204,0);stroke-width:2"); +// else +// mapper.map( +// poly, +// "fill-opacity:0.5;fill:rgb(230,104,58);stroke:rgb(230," +// "104,58);stroke-width:2"); +// } else { +// mapper.map(poly, +// "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153," +// "204,0);stroke-width:2"); +// } +// +// if (chain_d == 1) { +// auto lines = s_comp->get_level(1); +// for (int i = 0; i < lines.size(); ++i) { +// if (chain_v.coeffRef(i) != 0) { +// std::vector line_vec; +// for (auto k : lines.at(i)) { +// point_t pt = s_comp->get_points().at(k); +// point_type pt_g(pt[0], pt[1]); +// line_vec.push_back(pt_g); +// } +// polygon_type line_g(line_vec.begin(), line_vec.end()); +// mapper.map( +// line_g, +// "fill-opacity:0.5;fill:rgb(250,0,0);stroke:rgb(250," +// "0,0);stroke-width:2"); +// } +// } +// } +// } +// } chain_t bounding_chain::get_bounding_chain(chain_t& chain) { int chain_d; vector_t chain_v; std::tie(chain_d, chain_v) = chain; - if (chain_d >= this->s_comp->dimension()) throw NonZeroChain(); + if (chain_d >= s_comp->dimension()) throw non_zero_chain(); Eigen::LeastSquaresConjugateGradient lscg; - lscg.compute(*(this->boundary_matrices.at(chain_d))); + lscg.compute(*(boundary_matrices.at(chain_d))); vector_t bound_chain(round_vec(lscg.solve(chain_v))); + vector_t result(round_vec(*(boundary_matrices.at(chain_d)) * bound_chain)); - vector_t result( - round_vec(*(this->boundary_matrices.at(chain_d)) * bound_chain)); if (equals(result, chain_v)) return chain_t(chain_d - 1, bound_chain); else { - std::cout << "where to write the svg file?"; - std::string filename; - std::cin >> filename; - this->draw_chain(filename, chain_t(1, bound_chain)); - std::cout << '\n'; - std::cout << result << '\n'; - std::cout << chain_v << '\n'; - throw NonZeroChain(); + // std::cout << "where to write the svg file?"; + // std::string filename; + // std::cin >> filename; + // draw_chain(filename, chain_t(1, bound_chain)); + // std::cout << '\n'; + // std::cout << result << '\n'; + // std::cout << chain_v << '\n'; + throw non_zero_chain(); } } }; diff --git a/lib/scomplex/coeff_flow.cpp b/lib/scomplex/coeff_flow.hpp similarity index 66% rename from lib/scomplex/coeff_flow.cpp rename to lib/scomplex/coeff_flow.hpp index 44f084b..e4e613f 100644 --- a/lib/scomplex/coeff_flow.cpp +++ b/lib/scomplex/coeff_flow.hpp @@ -1,8 +1,8 @@ #include +#include #include -#include -#include #include +#include namespace gsimp { @@ -19,12 +19,13 @@ chain_t coeff_flow(simplicial_complex& s_comp, // if (chain_dim(p) != s_comp.dimension() - 1) throw out_of_context(); // 01 - chain_t c_chain(s_comp.dimension(), - vector_t(s_comp.get_level_size(s_comp.dimension()))); - std::vector seen_sigma(s_comp.get_level_size(s_comp.dimension()), - false); - std::vector seen_tau(s_comp.get_level_size(s_comp.dimension() - 1), - false); + chain_t c_chain = s_comp.new_chain(s_comp.dimension()); + + std::vector seen_sigma( + (size_t)s_comp.get_level_size(s_comp.dimension()), false); + + std::vector seen_tau( + (size_t)s_comp.get_level_size(s_comp.dimension() - 1), false); // 04 size_t tau_i, sigma_i(s_comp.cell_to_index(sigma_0)); @@ -38,7 +39,8 @@ chain_t coeff_flow(simplicial_complex& s_comp, // // 08 while (not queue.empty()) { - cell_t sigma, tau; + cell_t sigma; + cell_t tau; double c; std::tie(sigma, tau, c) = queue.front(); queue.pop(); @@ -48,7 +50,7 @@ chain_t coeff_flow(simplicial_complex& s_comp, // // 10 if (seen_sigma.at(sigma_i) && chain_val(c_chain, sigma_i) != c) throw no_bounding_chain(); - else if (not seen_sigma.at(sigma_i)) { + else if (not(seen_tau.at(tau_i) && seen_sigma.at(sigma_i))) { // 14 seen_tau.at(tau_i) = true; seen_sigma.at(sigma_i) = true; @@ -71,19 +73,32 @@ chain_t coeff_flow(simplicial_complex& s_comp, // int sigma_p_ind; cell_t sigma_p; std::tie(sigma_p_ind, sigma_p) = tau_cofaces[0]; // 20 - double c_p = - (sigma_p_ind * - (chain_val(p, tau_i) - - s_comp.boundary_inclusion_index(tau, sigma) * c)); - for (cell_t tau_p : s_comp.cell_boundary(sigma_p)) { - if (not seen_tau.at(s_comp.cell_to_index(tau_p))) { + double c_p = sigma_p_ind*(chain_val(p,tau_i)-s_comp.boundary_inclusion_index(tau, sigma)*c); + for (cell_t tau_p : s_comp.cell_boundary(sigma_p)) + if (not seen_tau.at(s_comp.cell_to_index(tau_p))) queue.emplace(sigma_p, tau_p, c_p); - seen_tau.at(tau_i) = true; - } - } } } } return c_chain; } + +chain_t coeff_flow_embedded(simplicial_complex& s_comp, chain_t p) { + if (chain_dim(p) != s_comp.dimension() - 1) throw out_of_context(); + + cell_t sigma; + double c; + auto level = s_comp.get_level(s_comp.dimension() - 1); + for (auto tau : s_comp.get_level(s_comp.dimension() - 1)) { + if (s_comp.get_cofaces(tau).size() < 2) { + int sigma_ind; + size_t tau_i = s_comp.cell_to_index(tau); + std::tie(sigma_ind, sigma) = s_comp.get_cof_and_ind(tau)[0]; + c = sigma_ind * chain_val(p, tau_i); + chain_t sol = coeff_flow(s_comp, p, sigma, c); + return sol; + } + } + throw out_of_context(); +} }; // namespace gsimp diff --git a/lib/scomplex/path_snapper.cpp b/lib/scomplex/path_snapper.cpp index 23ce003..7687e2a 100644 --- a/lib/scomplex/path_snapper.cpp +++ b/lib/scomplex/path_snapper.cpp @@ -106,7 +106,7 @@ chain_t path_snapper::snap_path_to_chain(std::vector path) { chain_t rep = p_impl->s_comp->new_chain(1); for (auto pair : pair_seq) { size_t index = p_impl->s_comp->cell_to_index(std::get<0>(pair)); - chain_val(rep, index) = std::get<1>(pair); + chain_val(rep, index) += std::get<1>(pair); } return rep; } @@ -127,10 +127,12 @@ std::vector path_snapper::point_sequence_to_index( } chain_t path_snapper::index_sequence_to_chain(std::vector ind_path) { - chain_t chain = p_impl->s_comp->new_chain(1); auto it = ind_path.begin(); - for (; std::next(it) != ind_path.end(); ++it) - chain_val(chain, p_impl->s_comp->cell_to_index({*it, *std::next(it)})); + chain_t chain = p_impl->s_comp->new_chain(1); + for (; std::next(it) != ind_path.end(); ++it) { + auto ind = p_impl->s_comp->cell_to_index({*it, *std::next(it)}); + chain_val(chain, ind) += (*it < *std::next(it)) ? 1 : -1; + } return chain; } diff --git a/lib/scomplex/simplicial_complex.cpp b/lib/scomplex/simplicial_complex.cpp index 14ce7c1..19ef127 100644 --- a/lib/scomplex/simplicial_complex.cpp +++ b/lib/scomplex/simplicial_complex.cpp @@ -157,6 +157,7 @@ std::vector> simplicial_complex::get_bdry_and_ind( }; + std::vector> simplicial_complex::get_bdry_and_ind_index( int d, size_t cell) { impl::simp_handle simp = p_impl->index_to_handle(d, cell); @@ -263,6 +264,7 @@ cell_t simplicial_complex::index_to_cell(int d, size_t ind) { auto sh = p_impl->levels.at(d)->at(ind); return p_impl->handle_to_cell(*sh); } + size_t simplicial_complex::cell_to_index(cell_t simp) { auto sh = p_impl->simplices.find(simp); return p_impl->simplices.key(sh); diff --git a/lib/scomplex/types.hpp b/lib/scomplex/types.hpp index 590267a..3db2d56 100644 --- a/lib/scomplex/types.hpp +++ b/lib/scomplex/types.hpp @@ -1,14 +1,14 @@ #pragma once +#include #include #include #include - namespace gsimp { // geometric types -typedef typename std::vector point_t; -typedef typename std::vector cell_t; +typedef std::vector point_t; +typedef std::vector cell_t; // linear algebra types typedef typename Eigen::SparseMatrix matrix_t; @@ -16,7 +16,48 @@ typedef typename Eigen::SparseVector vector_t; // chains typedef typename std::pair chain_t; + int& chain_dim(chain_t& p) { return std::get<0>(p); } vector_t& chain_rep(chain_t& p) { return std::get<1>(p); } double& chain_val(chain_t& p, size_t i) { return std::get<1>(p).coeffRef(i); } -}; // namespace gsimp + +chain_t add(chain_t chain1, chain_t chain2) { + if (std::get<0>(chain1) == std::get<0>(chain2)) + return chain_t(std::get<0>(chain1), + std::get<1>(chain1) + std::get<1>(chain2)); + std::cout << "chains must have the same dimension to be added (+)"; + throw std::exception(); +} + +void add_to(chain_t chain1, chain_t chain2) { + if (std::get<0>(chain1) == std::get<0>(chain2)) + std::get<1>(chain1) += std::get<1>(chain2); + else { + std::cout << "chains must have the same dimension to be added (+=)"; + throw std::exception(); + } +} + +chain_t subtract(chain_t chain1, chain_t chain2) { + if (std::get<0>(chain1) == std::get<0>(chain2)) + return chain_t(std::get<0>(chain1), + std::get<1>(chain1) - std::get<1>(chain2)); + std::cout << "chains must have the same dimension to be added (-)"; + throw std::exception(); +} + +void subtract_to(chain_t chain1, chain_t chain2) { + if (std::get<0>(chain1) == std::get<0>(chain2)) + std::get<1>(chain1) -= std::get<1>(chain2); + else { + std::cout << "chains must have the same dimension to be added (-=)"; + throw std::exception(); + } +} + +chain_t prod(double coef, chain_t chain) { + return chain_t(std::get<0>(chain), std::get<1>(chain) * coef); +} + +void prod_to(double coef, chain_t chain) { std::get<1>(chain) *= coef; } +} // namespace gsimp diff --git a/src/test0.cpp b/src/test0.cpp index 8637fd8..3dbf707 100644 --- a/src/test0.cpp +++ b/src/test0.cpp @@ -7,20 +7,21 @@ * 1-cell index, and outputs the cofaces of the given cell, (together with the * boundary index) and the faces of the cell. **/ +#include #include #include #include -#include #include +#include int my_char(gsimp::point_t point) { return 1; } int& do_this(int& a) { return a; } int main(int argc, char* argv[]) { - std::vector points_v; - std::vector cells_v; + std::vector points_v; + std::vector cells_v; std::string filename; if (argc == 1) { diff --git a/src/test1.cpp b/src/test1.cpp index 50c2431..960ec85 100644 --- a/src/test1.cpp +++ b/src/test1.cpp @@ -1,11 +1,20 @@ /** + * PATH SNAPPING TEST: + * takes as argument a path to a file containing the description of a + * triangulation as output by qhull (the last column of the point coordinates + * is ignored). + * comes up with a path in this complex (hard coded) and snaps it to the + * complex, the only point of this is to test the libraries for errors. **/ + #include #include #include #include #include + +// testing #include int my_char(gsimp::point_t point) { return 1; } diff --git a/src/test2.cpp b/src/test2.cpp new file mode 100644 index 0000000..e27c8df --- /dev/null +++ b/src/test2.cpp @@ -0,0 +1,60 @@ +/** + * PATH SNAPPING TEST: + * takes as argument a path to a file containing the description of a + * triangulation as output by qhull (the last column of the point coordinates + * is ignored). + * outputs the complex to an svg file called my_complex + **/ + +#include +#include + +#include +#include +#include +#include +#include + +int my_char(gsimp::point_t point) { return 1; } + +int main(int argc, char* argv[]) { + std::vector points_v; + std::vector cells_v; + + std::string filename; + if (argc == 1) { + std::cout << "insert file name: "; + std::cin >> filename; + } else + filename = argv[1]; + std::tie(points_v, cells_v) = parse_qhull_file(filename); + + gsimp::simplicial_complex s_comp(points_v, cells_v); + + gsimp::drawing drawing("my_complex", 500, 500); + + for (auto cell : s_comp.get_level(2)) { + std::vector> cell_pt; + for (auto i : cell) cell_pt.push_back(points_v[i]); + drawing.draw_triangle(cell_pt, "#00ff00"); + } + + for (auto cell : s_comp.get_level(1)) { + std::vector> cell_pt; + for (auto i : cell) cell_pt.push_back(points_v[i]); + drawing.draw_edge(cell_pt[0], cell_pt[1], "#000000", 0.002); + } + + gsimp::path_snapper p_snap(s_comp); + std::vector snapped = + p_snap.snap_path_to_indices({{-2, -2}, {-2, 2}, {2, 2}, {-2, -2}}); + auto it = snapped.begin(); + auto trail = it++; + for (; it != snapped.end(); ++it, ++trail) { + auto init = s_comp.get_points()[*trail]; + auto dest = s_comp.get_points()[*it]; + drawing.draw_edge(init, dest, "#ff0000", 0.003); + } + + return 0; +} diff --git a/src/test3.cpp b/src/test3.cpp new file mode 100644 index 0000000..e92babc --- /dev/null +++ b/src/test3.cpp @@ -0,0 +1,48 @@ +/** + * PATH SNAPPING TEST: + * takes as argument a path to a file containing the description of a + * triangulation as output by qhull (the last column of the point coordinates + * is ignored). + * outputs the complex to an svg file called my_complex + **/ + +#include +#include + +#include +#include +#include +#include +#include + +// testing +#include + +int my_char(gsimp::point_t point) { return 1; } + +int main(int argc, char* argv[]) { + std::vector points_v; + std::vector cells_v; + + std::string filename; + if (argc == 1) { + std::cout << "insert file name: "; + std::cin >> filename; + } else + filename = argv[1]; + std::tie(points_v, cells_v) = parse_qhull_file(filename); + + gsimp::simplicial_complex s_comp(points_v, cells_v); + + gsimp::drawing drawing("my_complex", 500, 500); + + gsimp::draw_complex(s_comp, drawing, 0.003, "#00ff00", "#000000"); + + { // scope for p_snap + gsimp::path_snapper p_snap(s_comp); + gsimp::draw_path(p_snap, drawing, {{-2, -2}, {2, -2}, {2, 2}, {-2, -2}}, + 0.003); + } + + return 0; +} diff --git a/src/test4.cpp b/src/test4.cpp new file mode 100644 index 0000000..6faebde --- /dev/null +++ b/src/test4.cpp @@ -0,0 +1,61 @@ +/** + * PATH SNAPPING TEST: + * takes as argument a path to a file containing the description of a + * triangulation as output by qhull (the last column of the point coordinates + * is ignored). + * outputs the complex to an svg file called my_complex + **/ + +#include +#include + +#include +#include +#include +#include +#include +#include + +// testing chain_calc +#include + +int my_char(gsimp::point_t point) { return 1; } + +int main(int argc, char* argv[]) { + std::vector points_v; + std::vector cells_v; + + std::string filename; + if (argc == 1) { + std::cout << "insert file name: "; + std::cin >> filename; + } else + filename = argv[1]; + std::tie(points_v, cells_v) = parse_qhull_file(filename); + + gsimp::simplicial_complex s_comp(points_v, cells_v); + gsimp::drawing drawing("my_complex", 500, 500); + gsimp::draw_complex(s_comp, drawing, 0.003, "#00ff00", "#000000"); + + std::vector path2 = {{-2, -2}, { 2, -2}, { 2, 2}, {-2, -2}}; + std::vector path = {{-2, -2}, {.3, .3}, {.2, -2}, {-2, -2}}; + + // auxiliary structures + gsimp::path_snapper p_snap(s_comp); + gsimp::bounding_chain c_calc(s_comp); + + gsimp::draw_path(p_snap, drawing, path, 0.003); + gsimp::draw_path(p_snap, drawing, path2, 0.003); + + gsimp::chain_t path_chain = p_snap.snap_path_to_chain(path); + gsimp::add_to(p_snap.snap_path_to_chain(path2),path_chain); + try { + gsimp::chain_t bound_chain = c_calc.get_bounding_chain(path_chain); + draw_2_chain(s_comp, drawing, bound_chain, + "#ff00ff' fill-opacity='0.8"); + } catch (const std::exception& e) { + std::cout << "oh well, nothing I can do"; + } + + return 0; +} diff --git a/src/test5.cpp b/src/test5.cpp new file mode 100644 index 0000000..3d1e72c --- /dev/null +++ b/src/test5.cpp @@ -0,0 +1,58 @@ +/** + * PATH SNAPPING TEST: + * takes as argument a path to a file containing the description of a + * triangulation as output by qhull (the last column of the point coordinates + * is ignored). + * outputs the complex to an svg file called my_complex + **/ + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +// testing chain_calc +#include + +int my_char(gsimp::point_t point) { return 1; } + +int main(int argc, char* argv[]) { + std::vector points_v; + std::vector cells_v; + + std::string filename; + if (argc == 1) { + std::cout << "insert file name: "; + std::cin >> filename; + } else + filename = argv[1]; + std::tie(points_v, cells_v) = parse_qhull_file(filename); + + gsimp::simplicial_complex s_comp(points_v, cells_v); + gsimp::drawing drawing("my_complex", 500, 500); + gsimp::draw_complex(s_comp, drawing, 0.003, "#00ff00", "#000000"); + + std::vector path2 = {{-2, -2}, { 2, -2}, { 2, 2}, {-2, -2}}; + std::vector path = {{-2, -2}, {.3, .3}, {.2, -2}, {-2, -2}}; + + // auxiliary structures + gsimp::path_snapper p_snap(s_comp); + gsimp::bounding_chain c_calc(s_comp); + + gsimp::draw_path(p_snap, drawing, path, 0.003); + gsimp::draw_path(p_snap, drawing, path2, 0.003); + + gsimp::chain_t path_chain = p_snap.snap_path_to_chain(path); + gsimp::add_to(p_snap.snap_path_to_chain(path2),path_chain); + + gsimp::chain_t bound_chain = gsimp::coeff_flow_embedded(s_comp,path_chain); + gsimp::draw_2_chain(s_comp, drawing, bound_chain, "#ff00ff' fill-opacity='0.8"); + + return 0; +} diff --git a/src/utils.hpp b/src/utils.hpp new file mode 100644 index 0000000..7348ac5 --- /dev/null +++ b/src/utils.hpp @@ -0,0 +1,2 @@ +#include +namespace gsimp