optimal_transport.h

Go to the documentation of this file.

/*
 *  Copyright (c) 2000-2022 Inria
 *  All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *
 *  * Redistributions of source code must retain the above copyright notice,
 *  this list of conditions and the following disclaimer.
 *  * Redistributions in binary form must reproduce the above copyright notice,
 *  this list of conditions and the following disclaimer in the documentation
 *  and/or other materials provided with the distribution.
 *  * Neither the name of the ALICE Project-Team nor the names of its
 *  contributors may be used to endorse or promote products derived from this
 *  software without specific prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 *
 *  Contact: Bruno Levy
 *
 *     https://www.inria.fr/fr/bruno-levy
 *
 *     Inria,
 *     Domaine de Voluceau,
 *     78150 Le Chesnay - Rocquencourt
 *     FRANCE
 *
 */
 
#ifndef H_EXPLORAGRAM_OPTIMAL_TRANSPORT_OPTIMAL_TRANSPORT_H
#define H_EXPLORAGRAM_OPTIMAL_TRANSPORT_OPTIMAL_TRANSPORT_H
 
#include <exploragram/basic/common.h>
 
struct NLMatrixStruct;
typedef NLMatrixStruct* NLMatrix;
 
namespace GEO {
    enum OTLinearSolver {
        OT_PRECG, OT_SUPERLU, OT_CHOLMOD
    };
}
 
#ifndef GOMGEN
 
#include <geogram/mesh/mesh.h>
#include <geogram/voronoi/RVD.h>
#include <geogram/delaunay/delaunay.h>
#include <geogram/NL/nl.h>
#include <geogram/NL/nl_matrix.h>
#include <geogram/third_party/HLBFGS/HLBFGS.h>
 
namespace GEO {
    class CentroidalVoronoiTesselation;
 
 
    class EXPLORAGRAM_API OptimalTransportMap {
    public:
    OptimalTransportMap(
        index_t dimension,
        Mesh* mesh,
        const std::string& delaunay = "default",
        bool BRIO = false
    );
 
    virtual ~OptimalTransportMap();
 
    index_t dimension() const {
        return dimension_;
    }
 
    Mesh& mesh() {
        return *mesh_;
    }
 
    void set_Newton(bool x) {
        newton_ = x;
    }
 
    void set_points(
        index_t nb_points, const double* points, index_t stride=0
    );
 
    void set_air_particles(
        index_t nb_air_particles, const double* air_particles, index_t stride,
        double air_fraction
    ) {
        nb_air_particles_ = nb_air_particles;
        air_particles_ = air_particles;
        air_particles_stride_ = (stride == 0) ? dimension_ : stride;
        air_fraction_ = air_fraction;
        // "continuous air" mode (air fraction declared without any air
        // particle).
        clip_by_balls_ = (nb_air_particles == 0) && (air_fraction != 0.0);
    }
 
    double air_fraction() const {
        return air_fraction_;
    }
 
    index_t nb_air_particles() const {
        return nb_air_particles_;
    }
 
    void set_nu(index_t i, double nu);
 
 
    void set_Laguerre_centroids(double* x) {
        Laguerre_centroids_ = x;
    }
 
    void set_epsilon(double eps) {
        epsilon_ = eps;
    }
 
 
    void set_linsolve_epsilon(double eps) {
        linsolve_epsilon_ = eps;
    }
 
    void set_linsolve_maxiter(index_t maxiter) {
        linsolve_maxiter_ = maxiter;
    }
 
    void set_linesearch_maxiter(index_t maxiter) {
        linesearch_maxiter_ = maxiter;
    }
 
    void set_linesearch_init_iter(index_t init_iter) {
        linesearch_init_iter_ = init_iter;
    }
 
    void set_regularization(double eps_reg) {
        epsilon_regularization_ = eps_reg;
    }
 
 
    void set_linear_solver(OTLinearSolver solver) {
        linear_solver_ = solver;
    }
 
    void optimize(index_t max_iterations);
 
    void set_verbose(bool x) {
        verbose_ = x;
    }
 
    void optimize_full_Newton(index_t max_iterations, index_t n=0);
 
    void optimize_level(index_t b, index_t e, index_t max_iterations);
 
    void optimize_levels(
        const vector<index_t>& levels, index_t max_iterations
    );
 
    index_t nb_points() const {
        return weights_.size();
    }
 
    const double* point_ptr(index_t i) const {
        geo_debug_assert(i < nb_points());
        return &(points_dimp1_[dimp1_ * i]);
    }
 
    double weight(index_t i) const {
        return weights_[i];
    }
 
    void set_weight(index_t i, double val) {
        weights_[i] = val;
    }
 
    double potential(index_t i) const {
        return points_dimp1_[dimp1_*i + dimension_];
    }
 
    static void funcgrad_CB(
        index_t n, double* x, double& f, double* g
    );
 
    static void newiteration_CB(
        index_t n, const double* x, double f, const double* g, double gnorm
    );
 
    RestrictedVoronoiDiagram* RVD() {
        return RVD_;
    }
 
    void set_save_RVD_iter(
        bool x,
        bool show_RVD_seed = false,
        bool last_iter_only = false
    ) {
        if(last_iter_only) {
            save_RVD_iter_ = false;
            save_RVD_last_iter_ = true;
        } else {
            save_RVD_iter_ = x;
        }
        show_RVD_seed_ = show_RVD_seed;
    }
 
    virtual void get_RVD(Mesh& M) = 0;
 
    virtual void compute_Laguerre_centroids(double* centroids) = 0;
 
    void update_sparsity_pattern();
 
    void new_linear_system(index_t n, double* x);
 
    void add_ij_coefficient(index_t i, index_t j, double a) {
        if(!user_H_g_) {
            nlAddIJCoefficient(i,j,a);
        } else {
            if(user_H_ != nullptr) {
                geo_debug_assert(user_H_->type == NL_MATRIX_SPARSE_DYNAMIC);
                nlSparseMatrixAdd((NLSparseMatrix*)user_H_, i, j, a);
            }
        }
    }
 
    void add_i_right_hand_side(index_t i, double a) {
        if(!user_H_g_) {
            nlAddIRightHandSide(i,a);
        }
    }
 
    void solve_linear_system();
 
    void set_initial_weight(index_t i, double w) {
        weights_[i] = w;
    }
 
    double total_mass() const {
        return total_mass_;
    }
 
    void compute_P1_Laplacian(
        const double* weights, NLMatrix Laplacian, double* measures
    );
 
    protected:
 
    double nu(index_t p) const {
        return nu_.size() == 0 ? constant_nu_ : nu_[p];
    }
 
 
    virtual void newiteration();
 
    void save_RVD(index_t id);
 
    void funcgrad(index_t n, double* w, double& f, double* g);
 
    virtual void call_callback_on_RVD() = 0;
 
    void eval_func_grad_Hessian(
        index_t n, const double* w,
        double& f, double* g
    );
 
    double gradient_threshold(index_t n) const {
        return ::sqrt(double(n) * geo_sqr(epsilon_ * constant_nu_));
    }
 
    public:
 
    class Callback {
    public:
        Callback(
            OptimalTransportMap* OTM
        ) : OTM_(OTM),
            Newton_step_(false),
            eval_F_(false),
            n_(0),
            w_(nullptr),
            g_(nullptr),
            mg_(nullptr) {
            weighted_ =
                OTM->mesh().vertices.attributes().is_defined("weight");
        }
 
        virtual ~Callback();
 
        void set_Laguerre_centroids(double* mg) {
            mg_ = mg;
        }
 
        bool has_Laguerre_centroids() const {
            return (mg_ != nullptr);
        }
 
        double* Laguerre_centroids() {
            return mg_;
        }
 
        void set_w(const double* w, index_t n) {
            w_ = w;
            n_ = n;
        }
 
        void set_Newton_step(bool Newton) {
            Newton_step_ = Newton;
        }
 
        bool is_Newton_step() const {
            return Newton_step_;
        }
 
        void set_nb_threads(index_t nb) {
            funcval_.assign(nb, 0.0);
        }
 
        void set_g(double* g) {
            g_ = g;
        }
 
        void set_eval_F(bool x) {
            eval_F_ = x;
        }
 
        double funcval() const {
            double result = 0.0;
            FOR(i,funcval_.size()) {
                result += funcval_[i];
            }
            return result;
        }
 
    protected:
        OptimalTransportMap* OTM_;
        bool weighted_;
        bool Newton_step_;
        bool eval_F_;
        vector<double> funcval_;
        index_t n_;
        const double* w_;
        double* g_;
        double* mg_;
    };
 
    protected:
    static OptimalTransportMap* instance_;
    index_t dimension_;
    index_t dimp1_; 
    Mesh* mesh_;
    Delaunay_var delaunay_;
    RestrictedVoronoiDiagram_var RVD_;
    vector<double> points_dimp1_;
    vector<double> weights_;
    double total_mass_;
    double constant_nu_; 
    vector<double> nu_;  
    double epsilon_;
    index_t current_call_iter_;
 
    Callback* callback_;
 
    std::string last_stats_;
    bool pretty_log_;
    index_t level_;
 
    bool save_RVD_iter_;
    bool save_RVD_last_iter_;
    bool show_RVD_seed_;
    index_t current_iter_;
    bool newton_;
    bool verbose_;
 
    double epsilon_regularization_;
 
    index_t nbZ_;
 
    double g_norm_;
 
    double measure_of_smallest_cell_;
 
    bool w_did_not_change_;
 
    double* Laguerre_centroids_;
 
    double linsolve_epsilon_;
 
    index_t linsolve_maxiter_;
 
    index_t linesearch_maxiter_;
 
    index_t linesearch_init_iter_;
 
    OTLinearSolver linear_solver_;
 
    const double* air_particles_;
 
    index_t nb_air_particles_;
 
    index_t air_particles_stride_;
 
    double air_fraction_;
 
    bool clip_by_balls_;
 
 
    bool user_H_g_;
 
    NLMatrix user_H_;
    };
 
}
 
#endif
 
#endif