/* This file is part of the Gudhi Library. The Gudhi library
* (Geometric Understanding in Higher Dimensions) is a generic C++
* library for computational topology.
*
* Author(s): Siargey Kachanovich
*
* Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France)
*
* 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 3 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, see .
*/
#ifndef GUDHI_WITNESS_COMPLEX_H_
#define GUDHI_WITNESS_COMPLEX_H_
#include
#include
#include
#include
#include "gudhi/reader_utils.h"
#include "gudhi/distance_functions.h"
#include "gudhi/Simplex_tree.h"
#include
#include
#include
#include
namespace Gudhi {
/*
template
class Simplex_tree;
*/
/*
typedef int Witness_id;
typedef int Landmark_id;
typedef int Simplex_handle; //index in vector complex_
*/
template
class Witness_complex: public Simplex_tree<> {
//class Witness_complex: public Simplex_tree {
//class Witness_complex {
private:
/*
template < typename Witness_id = int,
typename Landmark_id = int,
typename Simplex_handle = Simplex_handle >
struct Active_witness {
Witness_id witness_id;
Landmark_id landmark_id;
Simplex_handle simplex_handle;
};
*/
struct Active_witness {
int witness_id;
int landmark_id;
Simplex_handle simplex_handle;
Active_witness(int witness_id_, int landmark_id_, Simplex_handle simplex_handle_)
: witness_id(witness_id_),
landmark_id(landmark_id_),
simplex_handle(simplex_handle_)
{}
};
public:
//Simplex_tree<> st;
// typedef typename std::vector< Simplex_handle >::iterator Boundary_simplex_iterator;
// typedef boost::iterator_range Boundary_simplex_range;
// typedef typename std::vector< Simplex_handle >::iterator Skeleton_simplex_iterator;
// typedef boost::iterator_range< Skeleton_simplex_iterator > Skeleton_simplex_range;
// typedef IndexingTag Indexing_tag;
/** \brief Type for the value of the filtration function.
*
* Must be comparable with <. */
// typedef FiltrationValue Filtration_value;
/** \brief Key associated to each simplex.
*
* Must be a signed integer type. */
// typedef SimplexKey Simplex_key;
/** \brief Type for the vertex handle.
*
* Must be a signed integer type. It admits a total order <. */
typedef VertexHandle Vertex_handle;
/* Type of node in the simplex tree. */
typedef Simplex_tree_node_explicit_storage Node;
/* Type of dictionary Vertex_handle -> Node for traversing the simplex tree. */
typedef typename boost::container::flat_map Dictionary;
typedef typename Dictionary::iterator Simplex_handle;
/*
friend class Simplex_tree_node_explicit_storage< Simplex_tree >;
friend class Simplex_tree_siblings< Simplex_tree, Dictionary>;
friend class Simplex_tree_simplex_vertex_iterator< Simplex_tree >;
friend class Simplex_tree_boundary_simplex_iterator< Simplex_tree >;
friend class Simplex_tree_complex_simplex_iterator< Simplex_tree >;
friend class Simplex_tree_skeleton_simplex_iterator< Simplex_tree >;
*/
/* \brief Set of nodes sharing a same parent in the simplex tree. */
/* \brief Set of nodes sharing a same parent in the simplex tree. */
// typedef Simplex_tree_siblings Siblings;
typedef std::vector< double > Point_t;
typedef std::vector< Point_t > Point_Vector;
typedef std::vector< Vertex_handle > typeVectorVertex;
typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex;
typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool;
typedef int Witness_id;
typedef int Landmark_id;
typedef std::list< Vertex_handle > ActiveWitnessList;
/**
* /brief Iterative construction of the witness complex basing on a matrix of k nearest neighbours of the form {witnesses}x{landmarks}.
* Landmarks are supposed to be in [0,nbL-1]
*/
template< typename KNearestNeighbours >
void witness_complex(KNearestNeighbours & knn)
//void witness_complex(std::vector< std::vector< Vertex_handle > > & knn)
{
std::cout << "**Start the procedure witness_complex" << std::endl;
int k=2; /* current dimension in iterative construction */
//Construction of the active witness list
int nbW = knn.size();
int nbL = knn.at(0).size();
//VertexHandle vh;
typeVectorVertex vv;
typeSimplex simplex;
typePairSimplexBool returnValue;
int counter = 0;
/* The list of still useful witnesses
* it will diminuish in the course of iterations
*/
ActiveWitnessList active_w;// = new ActiveWitnessList();
for (int i=0; i != nbL; ++i) {
// initial fill of 0-dimensional simplices
// by doing it we don't assume that landmarks are necessarily witnesses themselves anymore
//vh = (Vertex_handle)i;
counter++;
vv = {i};
/* TODO Filtration */
//simplex = std::make_pair(vv, Filtration_value(0.0));
//returnValue = this->insert_simplex(simplex.first, simplex.second);
returnValue = insert_simplex(vv, Filtration_value(0.0));
/* TODO Error if not inserted : normally no need here though*/
}
//vv = {0};
//returnValue = insert_simplex(vv,Filtration_value(0.0));
std::cout << "Successfully added landmarks" << std::endl;
// PRINT2
print_sc(root()); std::cout << std::endl;
int u,v; // two extremities of an edge
if (nbL > 1) // if the supposed dimension of the complex is >0
/*
// THE BUGGY CODE
for (int i=0; i != nbW; ++i) {
// initial fill of active witnesses list
u = knn[i][0];
v = knn[i][1];
//Siblings * curr_sib = &root_;
//vh = (Vertex_handle)i;
vv = {u,v};
counter++;
returnValue = this->insert_simplex(vv,Filtration_value((double)counter));
//std::cout << "Null simplex is " << null_simplex()->first << std::endl;
if (returnValue.first != null_simplex())
{
active_w.push_back(*(new Active_witness(i,v,returnValue.first)));
}
for (typename ActiveWitnessList::iterator it1 = active_w.begin(); it1 != active_w.end(); ++it1)
std::cout << it1->simplex_handle->first << " ";
std::cout << std::endl;
//Simplex_handle sh = root_.members_.begin()+u;
//active_w.push_front(i);
}
*/
for (int i=0; i != nbW; ++i) {
// initial fill of active witnesses list
u = knn[i][0];
v = knn[i][1];
//Siblings * curr_sib = &root_;
//vh = (Vertex_handle)i;
vv = {u,v};
returnValue = this->insert_simplex(vv,Filtration_value(0.0));
print_sc(root()); std::cout << std::endl;
//std::cout << "Added edges" << std::endl;
}
//print_sc(root());
for (int i=0; i != nbW; ++i) {
// initial fill of active witnesses list
u = knn[i][0];
v = knn[i][1];
if ( u > v) {
u = v;
v = knn[i][0];
knn[i][0] = knn[i][1];
knn[i][1] = v;
}
Simplex_handle sh;
vv = {u,v};
sh = (root()->find(u))->second.children()->find(v);
active_w.push_back(i);
/*for (typename ActiveWitnessList::iterator it1 = active_w.begin(); it1 != active_w.end(); ++it1)
std::cout << it1->simplex->first << " ";
std::cout << std::endl; */
}
std::cout << "Successfully added edges" << std::endl;
while (!active_w.empty() && k+1 < nbL ) {
std::cout << "Started the step k=" << k << std::endl;
typename ActiveWitnessList::iterator it = active_w.begin();
while (it != active_w.end()) {
typeVectorVertex simplex_vector;
/* THE INSERTION: Checking if all the subfaces are in the simplex tree*/
// First sort the first k landmarks
int i = k;
int temp_swap;
VertexHandle inserted_vertex = knn[*it][k];
while (i>0 && knn[*it][i-1] > knn[*it][i])
{
temp_swap = knn[*it][i];
knn[*it][i] = knn[*it][i-1];
knn[*it][i-1] = temp_swap;
i--;
}
bool ok = all_faces_in(knn, *it, k, inserted_vertex);
if (ok)
{
for (i = 0; i != k+1; ++i)
{
simplex_vector.push_back(knn[*it][i]);
}
returnValue = insert_simplex(simplex_vector,0.0);
}
else
active_w.erase(it); //First increase the iterator and then erase the previous element
it++;
}
print_sc(root());
k++;
}
}
private:
void print_sc(Siblings * sibl)
{
if (sibl == NULL)
std::cout << "&";
else
print_children(sibl->members_);
}
void print_children(Dictionary map)
{
std::cout << "(";
if (!map.empty())
{
std::cout << map.begin()->first;
/*if (map.begin()->second.children() == root())
std::cout << "Sweet potato"; */
if (has_children(map.begin()))
print_sc(map.begin()->second.children());
typename Dictionary::iterator it;
for (it = map.begin()+1; it != map.end(); ++it)
{
std::cout << "," << it->first;
/*if (map.begin()->second.children() == root())
std::cout << "Sweet potato";*/
if (has_children(it))
print_sc(it->second.children());
}
}
std::cout << ")";
}
/** Check if all the facets of a simplex we are going to insert are in the simplex tree or not.
* The only purpose is to test if the witness is still active or not.
* Assuming here that the list of the first k witnessed landmarks is sorted
*/
template< typename KNearestNeighbours >
bool all_faces_in(KNearestNeighbours &knn, int witness_id, int k, VertexHandle inserted_vertex)
{
std::cout << "All face in with the landmark " << inserted_vertex << std::endl;
//std::list< VertexHandle > suffix;
Simplex_handle curr_sh = root()->find(knn[witness_id][0]);
Siblings * curr_sibl = root();
//VertexHandle curr_vh = curr_sh->first;
// CHECK ALL THE FACETS
Simplex_handle sh_bup = curr_sh; // the first vertex lexicographicly
for (int i = 0; i != k+1; ++i)
{
curr_sh = sh_bup;
if (knn[witness_id][i] != inserted_vertex)
{
for (int j = 0; j != k+1; ++j)
{
if (j != i)
{
if (curr_sibl->find(knn[witness_id][j]) == null_simplex())
return false;
curr_sh = curr_sibl->find(knn[witness_id][j]);
curr_sibl = self_siblings(curr_sh);
}//endif j!=i
}//endfor
}//endif
} //endfor
return true;
}
/**
* \brief Permutes the vector in such a way that the landmarks appear first
*/
void furthestPoints(Point_Vector &W, int nbP, std::string file_land, int dim, int nbL, Point_Vector &L)
{
//std::cout << "Enter furthestPoints "<< endl;
//Point_Vector *L = new Point_Vector();
double density = 5.;
int current_number_of_landmarks=0;
double curr_max_dist;
double curr_dist;
double mindist = 10005.;
int curr_max_w=0;
int curr_w=0;
srand(354698);
int rand_int = rand()% nbP;
//std::cout << rand_int << endl;
L.push_back(W[rand_int]);// first landmark is random
current_number_of_landmarks++;
while (1) {
curr_w = 0;
curr_max_dist = -1;
for(Point_Vector::iterator itW = W.begin(); itW != W.end(); itW++) {
//compute distance from w and L
mindist = 100000.;
for(Point_Vector::iterator itL = L.begin(); itL != L.end(); itL++) {
//curr_dist = distPoints(*itW,*itL);
curr_dist = euclidean_distance(*itW,*itL);
if(curr_dist < mindist) {
mindist = curr_dist;
}
}
if(mindist > curr_max_dist) {
curr_max_w = curr_w; //???
curr_max_dist = mindist;
}
curr_w++;
}
L.push_back(W[curr_max_w]);
current_number_of_landmarks++;
density = sqrt(curr_max_dist);
//std::cout << "[" << current_number_of_landmarks << ":" << density <<"] ";
if(L.size() == nbL) break;
}
//std::cout << endl;
return L;
}
}; //class Witness_complex
} // namespace Guhdi
#endif