skbl improve efficiency code constructor

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/trunk@467 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: c1feee162465b6ef0e00811c028d94177e4fa6fc

4 files changed, 154 insertions, 84 deletions

diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h
index 434538fe..10d64e98 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h
@@ -182,6 +182,15 @@ class Skeleton_blocker_simplex {
     return simplex_set.cend();
   }
 
+  typename std::set<T>::const_reverse_iterator rbegin() const {
+    return simplex_set.crbegin();
+  }
+
+  typename std::set<T>::const_reverse_iterator rend() const {
+    return simplex_set.crend();
+  }
+
+
   typename std::set<T>::iterator begin() {
     return simplex_set.begin();
   }
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Trie.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Trie.h
index 03650f73..f2a443dc 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Trie.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Trie.h
@@ -30,6 +30,8 @@
 
 #include <memory>
 #include <vector>
+#include <deque>
+#include <set>
 
 namespace Gudhi {
 
@@ -37,10 +39,10 @@ namespace Gudhi {
 namespace skbl {
 
 
-template<typename SkeletonBlockerComplex>
+template<typename SimplexHandle>
 struct Trie{
-	typedef typename SkeletonBlockerComplex::Vertex_handle Vertex_handle;
-	typedef typename SkeletonBlockerComplex::Simplex_handle Simplex_handle;
+	typedef SimplexHandle Simplex_handle;
+	typedef typename SimplexHandle::Vertex_handle Vertex_handle;
 
 	Vertex_handle v;
 	std::vector<std::shared_ptr<Trie> > childs;
@@ -65,6 +67,7 @@ public:
 			child->parent_ = this;
 		}
 	}
+
 	typedef typename Simplex_handle::Simplex_vertex_const_iterator Simplex_vertex_const_iterator;
 
 
@@ -192,6 +195,80 @@ public:
 		return stream;
 	}
 };
+
+
+template<typename SimplexHandle>
+struct Tries{
+	typedef typename SimplexHandle::Vertex_handle Vertex_handle;
+	typedef SimplexHandle Simplex_handle;
+
+	typedef Trie<Simplex_handle> STrie;
+
+
+	template<typename SimpleHandleOutputIterator>
+	Tries(unsigned num_vertices,SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end):
+	cofaces_(num_vertices,0){
+		for (auto i = 0u; i < num_vertices; ++i)
+			cofaces_[i] = new STrie(Vertex_handle(i));
+		for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) {
+			if (s_it->dimension() >= 1)
+				cofaces_[s_it->first_vertex()]->add_simplex(*s_it);
+		}
+	}
+
+	~Tries(){
+		for(STrie* t : cofaces_)
+			delete t;
+	}
+
+	//return a simplex that consists in all u such uv is an edge and u>v
+	Simplex_handle positive_neighbors(Vertex_handle v) const{
+		Simplex_handle res;
+		for(auto child : cofaces_[v]->childs)
+			res.add_vertex(child->v);
+		return res;
+	}
+
+	bool contains(const Simplex_handle& s) const{
+		auto first_v = s.first_vertex();
+		return cofaces_[first_v]->contains(s);
+	}
+
+	friend std::ostream& operator<<(std::ostream& stream, const Tries& tries){
+		for(auto trie : tries.cofaces_)
+			stream<<*trie<<std::endl;
+		return stream;
+	}
+
+	//init_next_dimension must be called first
+	std::vector<Simplex_handle> next_dimension_simplices() const{
+		std::vector<Simplex_handle> res;
+		while(!to_see_.empty() && to_see_.front()->simplex().dimension()==current_dimension_){
+			res.emplace_back(to_see_.front()->simplex());
+			for(auto child : to_see_.front()->childs)
+				to_see_.push_back(child.get());
+			to_see_.pop_front();
+		}
+		++current_dimension_;
+		return res;
+	}
+
+	void init_next_dimension() const{
+		for(auto trie : cofaces_)
+			to_see_.push_back(trie);
+	}
+
+private:
+	mutable std::deque<STrie*> to_see_;
+	mutable unsigned current_dimension_=0;
+
+
+	std::vector<STrie*> cofaces_;
+
+};
+
+
+
 }
 
 }
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h
index 1f5ca238..0b397f56 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h
@@ -68,7 +68,7 @@ class Simplex_around_vertex_iterator :
 	// Link_vertex_handle == Complex_Vertex_handle but this renaming helps avoiding confusion
 
 
-	typedef typename Gudhi::skbl::Trie<Complex> Trie;
+	typedef typename Gudhi::skbl::Trie<Simplex_handle> Trie;
 
 
 private:
@@ -76,7 +76,7 @@ private:
 	Vertex_handle v;
 	std::shared_ptr<Link> link_v;
 	std::shared_ptr<Trie> trie;
-	std::list<Trie*> nodes_to_be_seen; // todo regrouper
+	std::list<Trie*> nodes_to_be_seen; // todo deque
 
 public:
 	Simplex_around_vertex_iterator():complex(0){
@@ -199,6 +199,12 @@ public:
 		return first_node->simplex();
 	}
 
+//private:
+	Simplex_handle get_trie_address() const{
+		assert(!nodes_to_be_seen.empty());
+		return nodes_to_be_seen.front();
+	}
+
 private:
 	void set_end(){
 		nodes_to_be_seen.clear();
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
index 81fc28cb..893c65c6 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
@@ -173,23 +173,8 @@ public:
 	}
 
 private:
-	typedef Trie<Skeleton_blocker_complex<SkeletonBlockerDS>> STrie;
-
-	template<typename SimpleHandleOutputIterator>
-	/**
-	 * return a vector of simplex trie for every vertex
-	 */
-	std::vector<STrie*> compute_cofaces(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end) {
-		std::vector<STrie*> cofaces(num_vertices(), 0);
-		for (auto i = 0u; i < num_vertices(); ++i)
-			cofaces[i] = new STrie(Vertex_handle(i));
-		for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) {
-			if (s_it->dimension() >= 1)
-				cofaces[s_it->first_vertex()]->add_simplex(*s_it);
-		}
-		return cofaces;
-	}
-
+	//	typedef Trie<Skeleton_blocker_complex<SkeletonBlockerDS>> STrie;
+	typedef Trie<Simplex_handle> STrie;
 
 public:
 	/**
@@ -201,6 +186,16 @@ public:
 	: num_vertices_(0),
 	  num_blockers_(0),
 	  visitor(visitor_) {
+		add_vertex_and_edges(simplex_begin,simplex_end);
+
+		if (!is_flag_complex)
+			// need to compute blockers
+			add_blockers(simplex_begin,simplex_end);
+	}
+
+private:
+	template<typename SimpleHandleOutputIterator>
+	void add_vertex_and_edges(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end){
 		std::vector<std::pair<Vertex_handle, Vertex_handle>> edges;
 		// first pass to add vertices and edges
 		int num_vertex = -1;
@@ -212,71 +207,53 @@ public:
 
 		for (const auto& e : edges)
 			add_edge(e.first, e.second);
+	}
 
-		if (!is_flag_complex) {
-			// need to compute blockers
-			std::vector<STrie*> cofaces(compute_cofaces(simplex_begin, simplex_end));
-			std::vector<Simplex_handle> max_faces;
-
-			for (auto i = 0u; i < num_vertices(); ++i) {
-				auto max_faces_i = cofaces[i]->maximal_faces();
-				max_faces.insert(max_faces.end(), max_faces_i.begin(), max_faces_i.end());
-			}
-
-			// all time here
-
-			// for every maximal face s, it picks its first vertex v and check for all nv \in N(v)
-			// if s union nv is in the complex, if not it is a blocker.
-			for (auto max_face : max_faces) {
-				if (max_face.dimension() > 0) {
-					auto first_v = max_face.first_vertex();
-					for (auto nv : vertex_range(first_v)) {
-						if (!max_face.contains(nv) && max_face.first_vertex() < nv) {
-							// check that all edges in vertices(max_face)\cup nv are here
-							// since max_face is a simplex, we only need to check that edges
-							// (x nv) where x \in max_face are present
-							bool all_edges_here = true;
-							for (auto x : max_face)
-								if (!contains_edge(x, nv)) {
-									all_edges_here = false;
-									break;
-								}
-							if (all_edges_here) {  // eg  this->contains(max_face)
-								max_face.add_vertex(nv);
-								if (!cofaces[first_v]->contains(max_face)) {
-									// if there exists a blocker included in max_face, we remove it
-									// as it is not a minimum missing face
-									// the other alternative would be to check to all properfaces
-									// are in the complex before adding a blocker but that
-									// would be more expensive if there are few blockers
-									std::vector<Blocker_handle> blockers_to_remove;
-									bool is_blocker = true;
-									for (auto b : blocker_range(first_v))
-										if (b->contains(max_face)){
-											blockers_to_remove.push_back(b);
-										}
-										else
-											if (max_face.contains(*b)) {
-												is_blocker = false;
-												break;
-											}
-									for (auto b : blockers_to_remove)
-										this->delete_blocker(b);
-
-									if(is_blocker)
-										add_blocker(max_face);
-								}
-								max_face.remove_vertex(nv);
-							}
+	template<typename SimpleHandleOutputIterator>
+	void add_blockers(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end){
+		Tries<Simplex_handle> tries(num_vertices(),simplex_begin,simplex_end);
+		tries.init_next_dimension();
+		auto simplices(tries.next_dimension_simplices());
+
+		while(!simplices.empty()){
+			simplices = tries.next_dimension_simplices();
+			for(auto& sigma : simplices){
+				//common_positive_neighbors is the set of vertices u such that
+				//for all s in sigma, us is an edge and u>s
+				Simplex_handle common_positive_neighbors(tries.positive_neighbors(sigma.last_vertex()));
+				for(auto sigma_it = sigma.rbegin(); sigma_it != sigma.rend(); ++sigma_it/**/)
+					if(sigma_it != sigma.rbegin())
+						common_positive_neighbors.intersection(tries.positive_neighbors(*sigma_it));
+
+				for(auto x : common_positive_neighbors){
+					//first test that all edges sx are here for all s in sigma
+					bool all_edges_here = true;
+					for(auto s : sigma)
+						if(!contains_edge(x,s)){
+							all_edges_here = false;
+							break;
 						}
-					}
+					if(!all_edges_here) continue;
+
+					// all edges of sigma \cup x are here
+					//we have a blocker if all proper faces of sigma \cup x
+					// are in the complex and if sigma \cup x is not in the complex
+					//the first is equivalent at checking if blocks(sigma \cup x) is true
+					//as all blockers of lower dimension have already been computed
+					sigma.add_vertex(x);
+					if(!tries.contains(sigma)&&!blocks(sigma))
+						add_blocker(sigma);
+					sigma.remove_vertex(x);
 				}
+
 			}
-			for (auto i = 0u; i < num_vertices(); ++i)
-				delete(cofaces[i]);
 		}
 	}
 
+
+public:
+
+
 	// We cannot use the default copy constructor since we need
 	// to make a copy of each of the blockers
 
@@ -824,11 +801,12 @@ private:
 	 * is a face of sigma.
 	 */
 	bool blocks(const Simplex_handle & sigma) const {
-		for (auto blocker : const_blocker_range()) {
-			if (sigma.contains(*blocker))
-				return true;
-		}
+		for(auto s : sigma)
+			for (auto blocker : const_blocker_range(s))
+				if (sigma.contains(*blocker))
+					return true;
 		return false;
+
 	}
 
 	//@}