Skeleton_blocker_link_complex.h

 /*    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):       David Salinas
  *
  *    Copyright (C) 2014  INRIA Sophia Antipolis-Mediterranee (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 <http://www.gnu.org/licenses/>.
  */
#ifndef GUDHI_SKELETON_BLOCKERS_LINK_COMPLEX_H_
#define GUDHI_SKELETON_BLOCKERS_LINK_COMPLEX_H_

#include "gudhi/Utils.h"
#include "gudhi/Skeleton_blocker_complex.h"

namespace Gudhi{

namespace skbl {

template<class ComplexType> class Skeleton_blocker_sub_complex;


template<typename ComplexType>
class Skeleton_blocker_link_complex : public Skeleton_blocker_sub_complex<ComplexType>
{
    template <typename T> friend class Skeleton_blocker_link_superior;
    typedef typename ComplexType::Edge_handle Edge_handle;

    typedef typename ComplexType::boost_vertex_handle boost_vertex_handle;


private:

    bool only_superior_vertices_;

public:
    typedef typename ComplexType::Vertex_handle Vertex_handle;
    typedef typename ComplexType::Root_vertex_handle Root_vertex_handle;

    typedef typename ComplexType::Simplex_handle Simplex_handle;
    typedef typename ComplexType::Root_simplex_handle Root_simplex_handle;

    typedef typename ComplexType::Blocker_handle Blocker_handle;


    typedef typename ComplexType::Simplex_handle::Simplex_vertex_const_iterator Simplex_handle_iterator;
    typedef typename ComplexType::Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_handle_iterator;


    Skeleton_blocker_link_complex(bool only_superior_vertices=false):only_superior_vertices_(only_superior_vertices){
    }

    Skeleton_blocker_link_complex(const ComplexType & parent_complex,const Simplex_handle& alpha_parent_adress,bool only_superior_vertices = false)
    :only_superior_vertices_(only_superior_vertices) {
        if(!alpha_parent_adress.empty())
            build_link(parent_complex,alpha_parent_adress);

    }

    Skeleton_blocker_link_complex(const ComplexType & parent_complex, Vertex_handle a_parent_adress, bool only_superior_vertices = false)
    :only_superior_vertices_(only_superior_vertices){
        Simplex_handle alpha_simplex(a_parent_adress);
        build_link(parent_complex,alpha_simplex);
    }


    Skeleton_blocker_link_complex(const ComplexType & parent_complex, Edge_handle edge, bool only_superior_vertices = false)
    :only_superior_vertices_(only_superior_vertices){
        Simplex_handle alpha_simplex(parent_complex.first_vertex(edge),parent_complex.second_vertex(edge));
        build_link(parent_complex,alpha_simplex);
    }

protected:


    void compute_link_vertices(const ComplexType & parent_complex,const Simplex_handle& alpha_parent_adress,bool only_superior_vertices,bool is_alpha_blocker = false){
        if(alpha_parent_adress.dimension()==0)
            // for a vertex we know exactly the number of vertices of the link (and the size of the corresponding vector)
            // thus we call a specific function that will reserve a vector with appropriate size
            this->compute_link_vertices(parent_complex,alpha_parent_adress.first_vertex(),only_superior_vertices_);
        else{
            // we compute the intersection of neighbors of alpha and store it in link_vertices
            Simplex_handle link_vertices_parent;
            parent_complex.add_neighbours(alpha_parent_adress,link_vertices_parent,only_superior_vertices);
            // For all vertex 'v' in this intersection, we go through all its adjacent blockers.
            // If one blocker minus 'v' is included in alpha then the vertex is not in the link complex.
            for (auto v_parent : link_vertices_parent){
                bool new_vertex = true;
                for (auto beta : parent_complex.const_blocker_range(v_parent))
                {
                    if(!is_alpha_blocker || *beta!=alpha_parent_adress){
                        new_vertex = !(alpha_parent_adress.contains_difference(*beta,v_parent));
                        if(!new_vertex) break;
                    }
                }
                if (new_vertex)
                    this->add_vertex(parent_complex.get_id(v_parent));
            }
        }
    }


    void compute_link_vertices(const ComplexType & parent_complex, Vertex_handle alpha_parent_adress,bool only_superior_vertices){
        // for a vertex we know exactly the number of vertices of the link (and the size of the corresponding vector
        this->skeleton.m_vertices.reserve(parent_complex.degree(alpha_parent_adress));

        // For all vertex 'v' in this intersection, we go through all its adjacent blockers.
        // If one blocker minus 'v' is included in alpha then the vertex is not in the link complex.
        for (auto v_parent : parent_complex.vertex_range(alpha_parent_adress)){
            if(!only_superior_vertices || v_parent.vertex> alpha_parent_adress.vertex)
                this->add_vertex(parent_complex.get_id(v_parent));
        }

    }


    void compute_link_edges(const ComplexType & parent_complex,const Simplex_handle& alpha_parent_adress,bool is_alpha_blocker = false){
        Simplex_handle_iterator y_link, x_parent , y_parent;
        // ----------------------------
        // Compute edges in the link
        // -------------------------

        if(this->num_vertices()<=1) return;

        for (auto x_link = this->vertex_range().begin() ;
                x_link!= this->vertex_range().end();
                ++x_link
        )
        {
            for (auto y_link = x_link; ++y_link != this->vertex_range().end(); ){
                Vertex_handle x_parent = *parent_complex.get_address(this->get_id(*x_link));
                Vertex_handle y_parent = *parent_complex.get_address(this->get_id(*y_link));
                if (parent_complex.contains_edge(x_parent,y_parent) ){
                    // we check that there is no blocker subset of alpha passing trough x and y
                    bool new_edge=true;
                    for (auto blocker_parent : parent_complex.const_blocker_range(x_parent))
                    {
                        if(!is_alpha_blocker || *blocker_parent!=alpha_parent_adress){
                            if  (blocker_parent->contains(y_parent))
                            {
                                new_edge = !(alpha_parent_adress.contains_difference(*blocker_parent,x_parent,y_parent));
                                if (!new_edge) break;
                            }
                        }
                    }
                    if (new_edge)
                        this->add_edge(*x_link,*y_link);
                }
            }
        }
    }


    boost::optional<Vertex_handle> give_equivalent_vertex(const ComplexType & other_complex,
            Vertex_handle address) const{
        Root_vertex_handle id((*this)[address].get_id());
        return other_complex.get_address(id);
    }

    /*
     * compute the blockers of the link if is_alpha_blocker is false.
     * Otherwise, alpha is a blocker, and the link is computed in the complex where
     * the blocker is anticollapsed.
     */
    void compute_link_blockers(const ComplexType & parent_complex,const Simplex_handle& alpha_parent,bool is_alpha_blocker = false){

        for (auto x_link : this->vertex_range()){

            Vertex_handle x_parent = * this->give_equivalent_vertex(parent_complex,x_link);

            for (auto blocker_parent : parent_complex.const_blocker_range(x_parent)){
                if(!is_alpha_blocker || *blocker_parent!=alpha_parent){
                    Simplex_handle sigma_parent(*blocker_parent);

                    sigma_parent.difference(alpha_parent);

                    if (sigma_parent.dimension()>=2 && sigma_parent.first_vertex() == x_parent){
                        Root_simplex_handle sigma_id(parent_complex.get_id(sigma_parent));
                        auto  sigma_link = this->get_simplex_address(sigma_id);
                        // ie if the vertices of sigma are vertices of the link
                        if(sigma_link){
                            bool is_new_blocker = true;
                            for(auto a : alpha_parent){
                                for(auto eta_parent : parent_complex.const_blocker_range(a)){
                                    if(!is_alpha_blocker || *eta_parent!=alpha_parent){
                                        Simplex_handle eta_minus_alpha(*eta_parent);
                                        eta_minus_alpha.difference(alpha_parent);
                                        if(eta_minus_alpha != sigma_parent &&
                                                sigma_parent.contains_difference(*eta_parent,alpha_parent)){
                                            is_new_blocker = false;
                                            break;
                                        }
                                    }
                                }
                                if (!is_new_blocker) break;
                            }
                            if (is_new_blocker)
                                this->add_blocker(new Simplex_handle(*sigma_link));

                        }
                    }
                }
            }
        }
    }


public:
    void build_link(const ComplexType & parent_complex,const Simplex_handle& alpha_parent_adress,bool is_alpha_blocker = false)
    {
        assert(is_alpha_blocker || parent_complex.contains(alpha_parent_adress));
        compute_link_vertices(parent_complex,alpha_parent_adress,only_superior_vertices_);
        compute_link_edges(parent_complex,alpha_parent_adress,is_alpha_blocker);
        compute_link_blockers(parent_complex,alpha_parent_adress,is_alpha_blocker);
    }


    friend void build_link_of_blocker(
            const ComplexType & parent_complex,
            Simplex_handle& blocker,
            Skeleton_blocker_link_complex & result){
        assert(blocker.dimension()>=2);
        assert(parent_complex.contains_blocker(blocker));
        result.clear();
        result.build_link(parent_complex,blocker,true);
    }


};

}

}  // namespace GUDHI

#endif /* SKELETON_BLOCKERS_LINK_COMPLEX_H_ */