Paving.hpp

Aller à la documentation de ce fichier.

#ifndef PAVING_H
#define PAVING_H
#include "Paving_config.hpp"
#include <Grapic.h>
#include <utility>
#include <iostream>
#include <fstream>
#include <vector>
#include <cassert>
#include <cmath>
#include <algorithm>
#include <ctime>




template<int N>
class Point final : private std::vector<TYPECOORD> {
public:
    using std::vector<TYPECOORD>::begin;
    using std::vector<TYPECOORD>::end;
    using std::vector<TYPECOORD>::push_back;
    using std::vector<TYPECOORD>::size;
    using std::vector<TYPECOORD>::clear;
    using std::vector<TYPECOORD>::operator[];


    Point();

    Point(const Point<N>&);

    // Constructeur par transfert
    //Point(Point<N> &&);

    Point(const TYPEVAL&);

    Point(const std::vector<TYPECOORD>&);

    Point(const std::vector<TYPECOORD>&, const TYPEVAL&);

    ~Point();

    inline TYPEVAL value() const {
        return this->val;
    }

    inline void setValue(const TYPEVAL &v) {
        this->val = v;
    }


    static double getDistance(const Point<N>&, const Point<N>&);


    Point<N> getMiddle(const Point<N>&, const Point<N>&) const;


    static long double getDeterminant(const std::vector<Point<N>>&);

    friend std::ostream& operator<<(std::ostream& os, const Point<N>& p) {
        if(p.size()==0) {
            os << "Empty point (" << p.value() << ")";
        } else {
            os << "(";
            for(unsigned int i = 0; i < N; i++) {
                os << p.at(i);
                if(i != p.size()-1) os << ", ";
            }
            os << ")";
            if(p.value() != 0) os << " [" << p.value() << "]";
        }
        return os;
    }

    friend inline bool operator==(const Point<N>& x, const Point<N>& y) {
        for(int ii=0; ii<N; ii++) {
            if(x[ii]!=y[ii])
                return false;
        }
        return true;
    }

    friend inline bool operator!=(const Point<N>& x, const Point<N>& y) {
        return !(x==y);
    }

    friend inline Point<N> operator+(const Point<N>& x, const Point<N>& y) {
        Point<N> res;
        for(int ii=0; ii<N; ii++) {
            res.push_back(x[ii]+y[ii]);
        }
        return res;
    }

    friend inline Point<N> operator-(const Point<N>& x, const Point<N>& y) {
        Point<N> res;
        for(int ii=0; ii<N; ii++) {
            res.push_back(x[ii]-y[ii]);
        }
        return res;
    }


private:
    TYPEVAL val; 



    long double getDeterminant_rec(const std::vector<std::vector<TYPECOORD>>&) const;
};


// Constructeur par défaut
template<int N>
Point<N>::Point() : std::vector<TYPECOORD>(), val(0) {
}

// Constructeur par copie
template<int N>
Point<N>::Point(const Point<N> &p) : std::vector<TYPECOORD>() {
    if(p.size() != N) {
        std::cerr << "\033[31mERROR\033[0m : Copying Point : Point size (" << p.size() << ") does not correspond to template parameter Point<" << N << ">" << std::endl;
        return;
    }
    for(auto elm : p) {
        this->push_back(elm);
    }
    this->val = p.value();
}

// Constructeur par transfert
/*
template<int N>
Point<N>::Point(Point<N> &&p) {
        std::cerr << "\033[31mERROR\033[0m : Transfering Point : Point size (" << p.size() << ") does not correspond to template parameter Point<" << N << ">" << std::endl;

}*/

// Constructeur par valeur
template<int N>
Point<N>::Point(const TYPEVAL &v) : std::vector<TYPECOORD>(), val(v) {
}

// Constructeur par conversion de vector
template<int N>
Point<N>::Point(const std::vector<TYPECOORD> &v) : std::vector<TYPECOORD>(), val(0) {
    if(v.size() != N) {
        std::cerr << "\033[31mERROR\033[0m : Making Point : Vector size (" << v.size() << ") does not correspond to template parameter Point<" << N << ">" << std::endl;
        return;
    }
    for(unsigned int ii=0; ii < v.size(); ii++) {
        this->push_back(v[ii]);
    }
    this->val = 0;
}

// Constructeur par conversion de vector et valeur
template<int N>
Point<N>::Point(const std::vector<TYPECOORD> &v, const TYPEVAL &value) : std::vector<TYPECOORD>() {
    std::cout << "Construit Point vector dim " << N << std::endl;
    if(v.size() != N) {
        std::cerr << "\033[31mERROR\033[0m : Making Point : Vector size (" << v.size() << ") does not correspond to template parameter Point<" << N << ">" << std::endl;
        return;
    }
    for(auto ii : v) {
        this->push_back(ii);
    }
    this->val = value;
}

// Destructeur
template<int N>
Point<N>::~Point() {
}

// Calcule la distance entre deux points de même dimension
template<int N>
double Point<N>::getDistance(const Point<N> &a, const Point<N> &b) {
    double res=0;
    if(a.size() != b.size())
        return -1;
    for(unsigned int ii=0; ii<a.size(); ii++) {
        res+=(a[ii]-b[ii])*(a[ii]-b[ii]);
    }
    return sqrt(res);
}

// Calcule le milieu
template<int N>
Point<N> Point<N>::getMiddle(const Point<N> &a, const Point<N> &b) const {
    Point<N> res;
    if(a.size() != b.size())
        return res;
    for(unsigned int ii=0; ii<a.size(); ii++) {
        res.push_back((a[ii]+b[ii])/2);
    }
    return res;
}

// Calcule le determinant entre deux points de même dimension
template<int N>
long double Point<N>::getDeterminant(const std::vector<Point<N>> &pts) {
    std::vector<std::vector<TYPECOORD>> v_pts;
    for(auto p: pts) {
        v_pts.push_back(p);
    }
    return pts[0].getDeterminant_rec(v_pts);
}

// Methode récursive du calcul de determinants
template<int N>
long double Point<N>::getDeterminant_rec(const std::vector<std::vector<TYPECOORD>> &pts) const {
    if(pts.size() == 2) { // Si dimension 2 retourne le produit
        return pts[0][0]*pts[1][1]-pts[1][0]*pts[0][1];
    } else { // Sinon decoupe la matrice avec ses mineurs
        int sign=-1;
        double tmpRes = 0;
        for(unsigned int ii=0; ii<pts.size(); ii++) {
            sign*=-1; //Le signe change a chaque itération
            std::vector<std::vector<TYPECOORD>> matrix_minor(pts);
            matrix_minor.erase(matrix_minor.begin());
            for(unsigned int jj=0; jj<pts.size()-1; jj++) {
                matrix_minor[jj].erase(matrix_minor[jj].begin()+ii);
            }
            tmpRes += sign*pts[0][ii]*Point<N>::getDeterminant_rec(matrix_minor);
        }
        return tmpRes;
    }
}




template<int N>
class Simplex final : private std::vector<Point<N>> {

public:
    using std::vector<Point<N>>::begin;
    using std::vector<Point<N>>::end;
    using std::vector<Point<N>>::push_back;
    using std::vector<Point<N>>::size;
    using std::vector<Point<N>>::clear;
    using std::vector<Point<N>>::operator[];

    Simplex();


    Simplex(const std::vector<Point<N>>&);


    Simplex(const int &);


    Simplex(const int &, const std::vector<Point<N>>&);

    ~Simplex();


    inline int getIndice() const {
        return this->indice;
    }


    double getVolume() const;


    bool contain(const Point<N>&) const;


    TYPEVAL get(const Point<N>&) const;

    friend std::ostream& operator<<(std::ostream& os, const Simplex<N>& s) {
        if(s.size()==0) {
            os << "Empty simplex";
        } else {
            os << s.getIndice() << " {" << std::endl;
            for(auto pt: s) {
                os << pt << std::endl;
            }
            os << "}";
        }
        return os;
    }

private:
    int indice; 
};


// Constructeur par défaut
template<int N>
Simplex<N>::Simplex() : std::vector<Point<N>>(), indice(-1) {
}

// Constructeur par conversion
template<int N>
Simplex<N>::Simplex(const std::vector<Point<N>> &pts) : std::vector<Point<N>>(pts), indice(-1) {
}

// Constructeur par défaut
template<int N>
Simplex<N>::Simplex(const int &i) : std::vector<Point<N>>(), indice(i) {
}

// Constructeur par conversion
template<int N>
Simplex<N>::Simplex(const int &i, const std::vector<Point<N>> &pts) : std::vector<Point<N>>(pts), indice(i) {
}

// Destructeur
template<int N>
Simplex<N>::~Simplex() {
}

// Renvoie le volume du simplexe
template<int N>
double Simplex<N>::getVolume() const {
    // Calcule la matrice pour le determinant
    auto pointSimplexe = this->at(0);
    std::vector<Point<N>> matrix;
    for(unsigned int ind=1; ind<=N; ind++) {
        auto pointFace = this->at(ind);
        matrix.push_back(pointFace-pointSimplexe);
    }
    // Calcule le determinant
    auto res = Point<N>::getDeterminant(matrix);

    int fact=1;
    for(int ii=1; ii<=N; ii++)
        fact*=ii;
    // Calcule le resultat final
    res/=fact;
    // Renvoie la valeur absolue
    return res>0 ? res : -res;
}

// Vérifie si un point se trouve dans un simplexe donné
template<int N>
bool Simplex<N>::contain(const Point<N> &p) const {
    if(getVolume()==0) return false;
    for(auto pointSimplexe: *this) {
        std::vector<Point<N>> matrixFace, matrixSommet; // matrices à comparer après calcul de determinant
        // Crée les deux matrices
        for(auto pointFace: *this) {
            if(pointSimplexe == pointFace) continue;
            matrixFace.push_back(p-pointFace);
            matrixSommet.push_back(pointSimplexe-pointFace);
        }
        // Calcule les determinants
        auto detFace = p.getDeterminant(matrixFace), detSommet = p.getDeterminant(matrixSommet);
        // On compare les determinants
        if(detFace*detSommet < 0) {
            // S'ils  n'ont pas le même signe, le point n'est pas dans le simplexe
            return false;
        }
    }
    // S'ils ont tous le même signe, le point est dans le simplexe
    return true;
}

// Renvoie la valeur du point qui se trouve dans ce sommet
template<int N>
TYPEVAL Simplex<N>::get(const Point<N>& p) const {
    TYPEVAL res = 0;
    if(!contain(p)) {
        std::cerr << "\033[31m\033[31mERROR\033[0m\033[0m: GET : simplex " << *this << " does not contain point " << p << std::endl;
        return 0;
    }
    //calcul volume simplex Vd
    auto volumeD = this->getVolume();
    Simplex<N> simplexQ;
    for(auto pointSommet: *this) {
        for(auto pointFace: *this) {
            if(pointSommet == pointFace) continue;
            simplexQ.push_back(pointFace);
        }
        simplexQ.push_back(p);
        // calcule volume Vq
        auto volumeQ = simplexQ.getVolume();
        // calcule la coordonnée barycentrique au point courant
        auto lambda = volumeQ/volumeD;
        // ajoute au résultat la coordonnée barycentrique multipliée par la veleur au point
        res += lambda * pointSommet.value();
        simplexQ.clear();
    }
    return res;
}




template<int N>
class Paving final {

public:

    Paving();


    Paving(const int &, const int &, const int &);


    Paving(const char* &f);


    Paving(const std::vector<Point<N>> &);

    ~Paving();


    inline std::vector<Point<N>> getPoints() {
        return this->points;
    }
    inline std::vector<Point<N>> getTreatedPoints() {
        return this->treated_points;
    }
    inline std::vector<std::vector<int>> getSimplices() {
        return this->simplices;
    }


    bool init();


    std::vector<Simplex<N>> sliceHypercube(const std::vector<Point<N>>&hypercube, const Point<N>&mid);




    bool add(Point<N> &p);


    TYPEVAL get(const Point<N> &p);



    Simplex<N> getSimplex(const Point<N> &p) const;



    std::vector<Point<N>> getPointsFromSimplex(const std::vector<int>&) const;


    void showPoints() const;
    void showTreatedPoints() const;
    void showSimplices() const;

    void draw() const;

    void show() const;


private:
    std::vector<Point<N>> points; 
    std::vector<Point<N>> treated_points; 
    std::vector<std::vector<int>> simplices; 

    Point<N> p_min, p_max;
    std::vector<Point<N>> user_points; 


    int getNearestPoint(const Point<N>&, const std::vector<Point<N>>&) const;



};


// Constructeur par défaut
template<int N>
Paving<N>::Paving() {
}

template<int N>
Paving<N>::Paving(const std::vector<Point<N>> &v) : points(v) {
}

template<int N>
Paving<N>::Paving(const int & n, const int & minv, const int & maxv) : p_min(Point<N>()), p_max(Point<N>())  {
    std::cout << "Generating " << n << " random points of dimension " << N << std::endl;
    for(int ii = 0; ii < N; ii++) {
        p_min.push_back(9999999);
        p_max.push_back(-9999999);
    }

    srand(time(NULL));
    for(int ii = 0; ii < n; ii++) {
        Point<N> pt((rand()%maxv)+ii);
        for(int jj=0; jj < N; jj++) {
            TYPEVAL val = (rand()%(maxv-minv))+minv;
            pt.push_back(val);
            p_min[jj] = val < p_min[jj] ? val : p_min[jj];
            p_max[jj] = val > p_max[jj] ? val : p_max[jj];
        }
        this->points.push_back(pt);
        pt.clear();
    }
}

// Constructeur avec nom de fichier
template<int N>
Paving<N>::Paving(const char* & f) : p_min(Point<N>()), p_max(Point<N>()) {
    std::cout << "Current data dimension " << N << std::endl << std::endl;
    std::cout << "Building Paving from file '" << f << "', " ;
    std::ifstream inputFile(f);
    int nbPoints, dim;
    for(int ii = 0; ii < N; ii++) {
        p_min.push_back(9999999);
        p_max.push_back(-9999999);
    }

    if(inputFile.is_open()) {
        inputFile >> nbPoints;
        inputFile >> dim;
        std::cout << "data dimension " << dim  << std::endl;
        if(dim != N) {
            std::cerr << "\033[31mERROR\033[0m : Reading file : Data dimension (" << dim << ") does not correspond to template parameter Paving<" << N << ">" << std::endl;
            return;
        }
        for(int ii=1; ii<=nbPoints; ii++) {
            Point<N> point;
            TYPECOORD coord;
            TYPEVAL val;
            for(int jj=0; jj<N; jj++) {
                inputFile >> coord;
                point.push_back(coord);
                p_min[jj] = val < p_min[jj] ? val : p_min[jj];
                p_max[jj] = val > p_max[jj] ? val : p_max[jj];
            }
            inputFile >> val;
            point.setValue(val);
            this->points.push_back(point);
        }
        inputFile.close();
    } else {
        std::cerr << "\033[31mERROR\033[0m: failed open file" << std::endl;
    }
}

// Destructeur
template<int N>
Paving<N>::~Paving() {
}

// INIT VERSION (D+1)-SIMPLEXE
template<int N>
bool Paving<N>::init() {
    std::cout << std::endl << "Paving initialisation" << std::endl << std::endl;
    if(N==2) grapic::winInit("PavIn", WINDIM_X, WINDIM_Y);

    auto extrem_coord = std::max(abs(*std::max_element(p_max.begin(),p_max.end())), abs(*std::min_element(p_min.begin(),p_min.end())));
    auto coeff = extrem_coord/100;
    if(coeff < 1) coeff = 1;

    //CREATION DU (N+1)-SIMPLEXE ENGLOBANT
    std::cout << "Creating englobing simplex..." << std::endl << std::endl;
    std::vector<Point<N>> tmp_points(points);
    points = std::vector<Point<N>>();
    Simplex<N+1> eng_simplex;
    Point<N+1> p;
    for(int ii = 0; ii < N+2; ii++) {
        for(int jj = 0; jj < N+1; jj++) {
            if(ii == N+1) {
                p.push_back(-extrem_coord);
            } else {
                if(ii==jj) p.push_back(extrem_coord);
                else p.push_back(0);
            }
        }
        eng_simplex.push_back(p);
        p.clear();
    }

    std::cout << "Start englobing simplex : " << std::endl;
    std::cout << eng_simplex << std::endl << std::endl;

    while(!tmp_points.empty()) {
        // Agrandissement du (N+1)-Simplexe
        for(int ii = 0; ii < N+1; ii++) {
            eng_simplex[ii][ii]+=coeff;
        }
        for(int ii = 0; ii < N+1; ii++) {
            eng_simplex[N+1][ii]-=coeff;
        }

        // Ajout des points restants à ce (N+1)-Simplexe
        for(unsigned int ii = 0; ii < tmp_points.size(); ii++) {
            Point<N> pt = tmp_points[ii];
            pt.push_back(0);
            Point<N+1> point;
            for(auto c: pt) point.push_back(c);
            point.push_back(0);
            if(eng_simplex.contain(point)) {
                //std::cout << point << " in " << eng_simplex << std::endl;
                points.push_back(tmp_points[ii]);
                tmp_points.erase(tmp_points.begin()+ii);
            }
        }
    }
    // Projection du (N+1)-Simplexe en dimension N
    Simplex<N> main_simplex;
    for(unsigned int ii = 0; ii < eng_simplex.size(); ii++) {
        if(ii == eng_simplex.size()-2) continue;
        Point<N+1> point = eng_simplex[ii];
        Point<N> pt;
        for(int jj = 0; jj < N; jj++) {
            pt.push_back(point[jj]);
        }
        main_simplex.push_back(pt);
    }


    std::cout << "Pre-refine englobing simplex : " << std::endl;
    std::cout << main_simplex << std::endl << std::endl;

    // AFFINAGE DU N-SIMPLEXE
    std::cout << "Refining englobing simplex..." << std::endl << std::endl;

    bool res = true;
    // Affinage du N-Simplexe en ramenant les points (t, 0, 0 ...), (0, t, 0 ...) ...
    while(res) {
        // Diminution du N-Simplexe
        for(int ii = 0; ii < N; ii++)
            main_simplex[ii][ii]-=coeff;
        // Teste si tous les points sont toujours dans le simplexe
        for(auto pt: points)
            if(!main_simplex.contain(pt))
                res = false;
    }
    // Remise en l'état
    for(int ii = 0; ii < N; ii++)
        main_simplex[ii][ii]+=coeff;
    p_max.clear();
    p_max.push_back(main_simplex[0][0]);
    p_max.push_back(main_simplex[0][0]);

    res = true;
    // Affinage du N-Simplexe en ramenant le point (-t, -t, -t ...)
    while(res) {
        for(int ii = 0; ii < N; ii++)
            main_simplex[N][ii]+=coeff;
        for(auto pt: points)
            if(!main_simplex.contain(pt))
                res = false;
    }
    for(int ii = 0; ii < N; ii++)
        main_simplex[N][ii]-=coeff;
    p_min.clear();
    auto v = main_simplex[N][0] > 0 ? 0 : main_simplex[N][0];
    p_min.push_back(v);
    p_min.push_back(v);

    std::cout << "Englobing simplex : " << std::endl;
    std::cout << main_simplex << std::endl << std::endl;

    show();

    // AJOUT DU SIMPLEXE ET DES POINTS
    std::vector<int> sim;
    for(unsigned int ii = 0; ii < main_simplex.size(); ii++) {
        treated_points.push_back(main_simplex[ii]);
        sim.push_back(ii);
    }
    simplices.push_back(sim);
    show();

    // Ajout de tous les points à l'intérieur de la boite
    Point<N> pt(main_simplex[main_simplex.size()-1]);
    res = true;
    while(!this->points.empty()) {
        int id = getNearestPoint(pt, this->points);
        pt = this->points[id];
        this->points.erase(this->points.begin()+id);
        if(!this->add(pt)) res = false;
        show();
    }
    user_points = std::vector<Point<N>>();
    std::cout << std::endl << std::endl << std::endl;
    return res;
}

// INIT VERSION HYPERCUBE
/*
template<int N>
void Paving<N>::init() {
    std::cout << std::endl << "Init Paving ..." << std::endl;
    grapic::winInit("PavIn", WINDIM_X, WINDIM_Y);

    // Création d'une boite englobante
    // Création du tableau des extremums, pour chaque axe de dimension, cherche le min et le max
    std::vector<Point<2>> extremums;
    for(int ii = 0; ii < N; ii++) {
        Point<2> extremum;
        extremum.push_back(999999);
        extremum.push_back(-999999);
        extremums.push_back(extremum);
    }
    for(auto pt: this->points) {
        for(unsigned int ii = 0; ii < pt.size(); ii++) {
            extremums[ii][0] = pt[ii] < extremums[ii][0] ? pt[ii] : extremums[ii][0];
            extremums[ii][1] = pt[ii] > extremums[ii][1] ? pt[ii] : extremums[ii][1];
        }
    }
    // Créeation des deux points extremaux, diagonale, et sert pour calculer un affichage correct
    for(auto extr: extremums) {
        p_min.push_back(extr[0]);
        p_max.push_back(extr[1]);
    }

    show();

    // Crée les coordonnées des points extremaux dans un vector a remettre dans l'ordre
    // Contient les valeurs issues de l'arbre binaire des possibilités des points extremaux
    bool b = true;
    std::vector<TYPECOORD> v_eng_points;
    for(unsigned int ii = 0; ii < extremums.size(); ii++) {
        int cpt = pow(2,N)/(2*(ii+1));
        for(int jj = 0; jj <= pow(2,N)-1; jj++) {
            if(jj%cpt == 0) b = b ? false : true;
            if(b) {
                v_eng_points.push_back(extremums[ii][1]);
            } else {
                v_eng_points.push_back(extremums[ii][0]);
            }
        }
    }

    // Crée les points extremaux à partir de l'arbre binaire, leur donne la valeur du point le plus proche
    std::vector<Point<N>> eng_points;
    int threshold = v_eng_points.size()/N;
    for(int ii = 0; ii < threshold; ii++) {
        Point<N> pt;
        for(int jj = 0; jj < N; jj++) {
            pt.push_back(v_eng_points[ii+jj*threshold]);
        }
        eng_points.push_back(pt);
    }
    std::cout << "Englobing points:" << std::endl;
    for(unsigned int ii = 0; ii < eng_points.size(); ii++) {
        eng_points[ii].setValue(points[getNearestPoint(eng_points[ii], this->points)].value());
        std::cout << eng_points[ii] << std::endl;
    }
    std::cout << std::endl;

    // Recherche du milieu, donne au milieu la valeur du point le plus proche
    Point<N> mid(p_min.getMiddle(p_min, p_max));
    mid.setValue(points[getNearestPoint(mid, this->points)].value());
    //std::cout << "Middle : " << mid << std::endl;

    // Création du contour de la boite (suite de points, pour éviter les simplexes qui se croisent), ajout au Pavage
    std::vector<Point<N>> cont_points;
    std::cout <<std::endl << "Outline:" << std::endl;
    cont_points.push_back(eng_points[0]);
    eng_points.erase(eng_points.begin());
    for(int ii = 0; !eng_points.empty(); ii++) {
        ii%=eng_points.size();
        for(int jj = 0; jj < N; jj++) {
            if(eng_points[ii][jj] == eng_points[eng_points.size()-1][jj]) {
                //Deux points sur le meme axe
                cont_points.push_back(eng_points[ii]);
                eng_points.erase(eng_points.begin()+ii);
                break;
            }
        }
    }
    for(auto pt: cont_points) {
        std::cout << pt << " -> ";
    }
    std::cout << std::endl << std::endl;
    this->treated_points = std::vector<Point<N>>(cont_points);
    this->treated_points.push_back(mid);

    // Création des nouveaux simplexes par combinaison de N points à la suite du contour, et du milieu
    std::vector<Simplex<N>> simps;
    for(unsigned int ii = 0; ii < cont_points.size(); ii++) {
        std::vector<int> simp;
        for(int jj = 0; jj < N; jj++) {
            simp.push_back((ii+jj)%cont_points.size());
        }
        simp.push_back(cont_points.size());

        simplices.push_back(simp);
    }
    showSimplices();
     show();

    // Ajout de tous les points à l'intérieur de la boite
    Point<N> p(mid);
     while(!this->points.empty()) {
         int id = getNearestPoint(p, this->points);
         p = this->points[id];
         this->points.erase(this->points.begin()+id);
         this->add(p);
         show();
     }
}
*/

// ADD
template<int N>
bool Paving<N>::add(Point<N> &p) {
    std::cout << "Adding point " << p << ": ";
    // Cherche le simplexe contenant
    Simplex<N> simplex = getSimplex(p);

    if(simplex.size()!=0) {
        // Ajoute le point aux points traités
        auto newInd = this->treated_points.size();
        this->treated_points.push_back(p);
        for(auto indSimplex: this->simplices[simplex.getIndice()]) {
            std::vector<int> sim;
            sim.push_back(newInd);
            for(auto indFace: this->simplices[simplex.getIndice()]) {
                if(indSimplex == indFace) continue;
                sim.push_back(indFace);
            }
            this->simplices.push_back(sim);
        }
        this->simplices.erase(this->simplices.begin()+simplex.getIndice());
        std::cout << "\033[32madded\033[0m" << std::endl;
        return true;
    } else {
        std::cout << "\033[31mnot in any simplex\033[0m" << std::endl;
        return false;
    }

}

// GET
template<int N>
TYPEVAL Paving<N>::get(const Point<N> &p) {
    TYPEVAL res = 0;
    std::cout << "Getting " << p << std::endl << std::endl;
    // Cherche le simplexe contenant
    Simplex<N> simplexD = getSimplex(p);
    std::cout << "Containing simplex: " << simplexD << std::endl;
    // Si un simplexe est trouvé
    if(simplexD.size()!=0) {
        res = simplexD.get(p);
        user_points.push_back(p);
    } else {
        std::cout << "\033[31mPoint " << p << " not in any simplex\033[0m" << std::endl;
    }
    return res;
}

// Point le plus proche
template<int N>
int Paving<N>::getNearestPoint(const Point<N> &p, const std::vector<Point<N>> &v) const {
    int res, ii=0;
    double d = 9999999, tmpD;
    for(auto elm : v) {
        tmpD = Point<N>::getDistance(p, elm);
        if(tmpD < d && p!=elm) {
            d = tmpD;
            res = ii;
        }
        ii++;
    }
    return res;
}

// Renvoie le simplexe contenant le point
template<int N>
Simplex<N> Paving<N>::getSimplex(const Point<N> &p) const {
    int ind=0;
    for(auto vec: this->simplices) {
        Simplex<N> simp(ind, getPointsFromSimplex(vec));
        if(simp.getVolume() > 0) {
            if(simp.contain(p)) {
                return simp;
            }
        }
        ind++;
    }
    return Simplex<N>();
}

// Renvoie les points d'un simplexe donné
template<int N>
std::vector<Point<N>> Paving<N>::getPointsFromSimplex(const std::vector<int> &s) const {
    std::vector<Point<N>> res;
    for(auto ind: s) {
        Point<N> pt = this->treated_points[ind];
        res.push_back(pt);
    }
    return res;
}

// showPoints
template<int N>
void Paving<N>::showPoints() const {
    std::cout << "Points (" << this->points.size() << ") :" << std::endl;
    for(auto p : this->points) {
        std::cout << "    " << p << std::endl;
    }
    std::cout << std::endl;
}

// showTreatedPoints
template<int N>
void Paving<N>::showTreatedPoints() const {
    std::cout << "Treated points (" << this->treated_points.size() << ") :" <<  std::endl;
    for(auto p : this->treated_points) {
        std::cout <<  "    " << p << std::endl;
    }
    std::cout << std::endl;
}

// showSimplices
template<int N>
void Paving<N>::showSimplices() const {
    std::cout << "Simplices (" << this->simplices.size() << ") :" << std::endl;
    int ind = 0;
    for(auto v: this->simplices) {
        Simplex<N> sim(ind++, getPointsFromSimplex(v));
        std::cout << sim << std::endl;
    }
    std::cout << std::endl;
}

// Dessin
template<int N>
void Paving<N>::draw() const {
    TYPECOORD height = p_max[0]-p_min[0];
    TYPECOORD width = p_max[1]-p_min[1];
    double coeff = WINDIM_X/width < WINDIM_Y/height ? (WINDIM_X)/width : (WINDIM_Y)/height;
    grapic::color(242,224,26);
    for(auto pt: this->points) {
        grapic::point((pt[0]-p_min[0])*coeff, (pt[1]-p_min[1])*coeff);
    }

    grapic::color(74, 239, 231);
    if(simplices.empty()) return;
    for(auto vec: this->simplices) {
        Point<N> a(this->treated_points[vec[0]]), b(this->treated_points[vec[1]]), c(this->treated_points[vec[2]]);
        grapic::line((a[0]-p_min[0])*coeff, (a[1]-p_min[1])*coeff, (b[0]-p_min[0])*coeff, (b[1]-p_min[1])*coeff);
        grapic::line((a[0]-p_min[0])*coeff, (a[1]-p_min[1])*coeff, (c[0]-p_min[0])*coeff, (c[1]-p_min[1])*coeff);
        grapic::line((c[0]-p_min[0])*coeff, (c[1]-p_min[1])*coeff, (b[0]-p_min[0])*coeff, (b[1]-p_min[1])*coeff);
    }

    grapic::color(239, 74, 1);
    if(user_points.empty()) return;
    Point<N> pt(user_points[user_points.size()-1]);
    Simplex<N> sim(getSimplex(pt));
    for(auto pt_sim: sim) {
        grapic::line((pt[0]-p_min[0])*coeff, (pt[1]-p_min[1])*coeff, (pt_sim[0]-p_min[0])*coeff, (pt_sim[1]-p_min[1])*coeff);
    }
}

// Affichage
template<int N>
void Paving<N>::show() const {
    if(N==2) {
        grapic::backgroundColor(0,0,0);
        grapic::winClear();
        draw();
        grapic::winDisplay();
        std::cout << "Press ENTER to continue" << std::endl;
        std::cin.ignore();
    }
}

#endif