GCC Code Coverage Report

Line	Branch	Exec	Source
1			/*
2			* Copyright (c) 2000-2022 Inria
3			* All rights reserved.
4			*
5			* Redistribution and use in source and binary forms, with or without
6			* modification, are permitted provided that the following conditions are met:
7			*
8			* * Redistributions of source code must retain the above copyright notice,
9			* this list of conditions and the following disclaimer.
10			* * Redistributions in binary form must reproduce the above copyright notice,
11			* this list of conditions and the following disclaimer in the documentation
12			* and/or other materials provided with the distribution.
13			* * Neither the name of the ALICE Project-Team nor the names of its
14			* contributors may be used to endorse or promote products derived from this
15			* software without specific prior written permission.
16			*
17			* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18			* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19			* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20			* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21			* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22			* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23			* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24			* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25			* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26			* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27			* POSSIBILITY OF SUCH DAMAGE.
28			*
29			* Contact: Bruno Levy
30			*
31			* https://www.inria.fr/fr/bruno-levy
32			*
33			* Inria,
34			* Domaine de Voluceau,
35			* 78150 Le Chesnay - Rocquencourt
36			* FRANCE
37			*
38			*/
39
40			#ifndef H_EXPLORAGRAM_OPTIMAL_TRANSPORT_OPTIMAL_TRANSPORT_H
41			#define H_EXPLORAGRAM_OPTIMAL_TRANSPORT_OPTIMAL_TRANSPORT_H
42
43			#include <exploragram/basic/common.h>
44
45			struct NLMatrixStruct;
46			typedef NLMatrixStruct* NLMatrix;
47
48			namespace GEO {
49			/**
50			* \brief Specifies the linear solver to be used
51			* with OptimalTransport.
52			*/
53			enum OTLinearSolver {
54			OT_PRECG, OT_SUPERLU, OT_CHOLMOD
55			};
56			}
57
58			#ifndef GOMGEN
59
60			#include <geogram/mesh/mesh.h>
61			#include <geogram/voronoi/RVD.h>
62			#include <geogram/delaunay/delaunay.h>
63			#include <geogram/NL/nl.h>
64			#include <geogram/NL/nl_matrix.h>
65			#include <geogram/third_party/HLBFGS/HLBFGS.h>
66
67			/**
68			* \file exploragram/optimal_transport/optimal_transport.h
69			* \brief Base class for semi-discrete optimal transport.
70			*/
71
72			namespace GEO {
73			class CentroidalVoronoiTesselation;
74
75
76			/**
77			* \brief Computes semi-discrete optimal transport maps.
78			* \details Computes an optimal transport map between two
79			* distributions. The first distribution is represented
80			* by a simplicial mesh. The second distribution is a sum
81			* of Diracs.
82			* The algorithm is described in the following references:
83			* - 3D algorithm: http://arxiv.org/abs/1409.1279
84			* - Earlier 2D version by Quentin M\'erigot:
85			* Q. Merigot. A multiscale approach to optimal transport.
86			* Computer Graphics Forum 30 (5) 1583--1592, 2011 (Proc SGP 2011).
87			* - Earlier article on OT and power diagrams:
88			* F. Aurenhammer, F. Hoffmann, and B. Aronov. Minkowski-type theorems
89			* and least-squares clustering. Algorithmica, 20:61-76, 1998.
90			*/
91			class EXPLORAGRAM_API OptimalTransportMap {
92			public:
93			/**
94			* \brief OptimalTransportMap constructor.
95			* \param[in] mesh the source distribution, represented as a nd mesh.
96			* \param[in] delaunay factory name of the Delaunay triangulation.
97			* \param[in] BRIO true if vertices are already ordered using BRIO
98			*/
99			OptimalTransportMap(
100			index_t dimension,
101			Mesh* mesh,
102			const std::string& delaunay = "default",
103			bool BRIO = false
104			);
105
106			/**
107			* \brief OptimalTransportMap destructor.
108			*/
109			virtual ~OptimalTransportMap();
110
111			/**
112			* \brief Gets the dimension.
113			* \return 2 for 2d, 3 for 3d.
114			*/
115			index_t dimension() const {
116		✗	return dimension_;
117			}
118
119			/**
120			* \brief Gets the mesh.
121			* \return a reference to the mesh
122			*/
123			Mesh& mesh() {
124		✗	return *mesh_;
125			}
126
127			/**
128			* \brief Sets whether Newton algorithm should be used.
129			* \details It is (for now) incompatible with multilevel.
130			* \param[in] x if set, Newton algorithm is used instead
131			* of BFGS.
132			*/
133			void set_Newton(bool x) {
134		✗	newton_ = x;
135			}
136
137			/**
138			* \brief Sets the points that define the target distribution.
139			* \details If air particles are used, then set_air_particles() needs
140			* to be called before set_points().
141			* \param[in] nb_points number of points in the target distribution
142			* \param[in] points an array of size nb_points * dimension() with the
143			* coordinates of the Diracs centers in the target
144			* distribution.
145			* \param[in] stride number of doubles between two consecutive points.
146			* If 0 (default), then point coordinates are considered to be packed.
147			*/
148			void set_points(
149			index_t nb_points, const double* points, index_t stride=0
150			);
151
152			/**
153			* \brief Sets the air particles that define the volume occupied by the
154			* free space.
155			* \details If air particles are used, then set_air_particles() needs
156			* to be called before set_points().
157			* \param[in] nb_air_particles number of air particles.
158			* \param[in] air_particles a pointer to the array of doubles with the
159			* coordinates of the air particles.
160			* \param[in] stride number of doubles between two consecutive air
161			* particles in the array, or 0 if tightly packed.
162			* \param[in] air_fraction the fraction of the total mass occupied by air.
163			*/
164			void set_air_particles(
165			index_t nb_air_particles, const double* air_particles, index_t stride,
166			double air_fraction
167			) {
168		✗	nb_air_particles_ = nb_air_particles;
169		✗	air_particles_ = air_particles;
170		✗	air_particles_stride_ = (stride == 0) ? dimension_ : stride;
171		✗	air_fraction_ = air_fraction;
172			// "continuous air" mode (air fraction declared without any air
173			// particle).
174		✗	clip_by_balls_ = (nb_air_particles == 0) && (air_fraction != 0.0);
175			}
176
177			/**
178			* \brief Gets the air fraction.
179			* \return the air fraction previously specified by set_air_particles()
180			*/
181			double air_fraction() const {
182		✗	return air_fraction_;
183			}
184
185			/**
186			* \brief Gets the number of air particles.
187			* \return the air fraction previously specified by set_air_particles()
188			*/
189			index_t nb_air_particles() const {
190		✗	return nb_air_particles_;
191			}
192
193			/**
194			* \brief Sets the desired mass at one of the Diracs.
195			* \details If unspecified, then default value is total mass
196			* divided by number of point. Note that the sum of all specified
197			* masses should match the total mass.
198			* \param[in] i the index of the dirac, in 0 .. nb_points-1
199			* \param[in] nu the desired mass at point i
200			*/
201			void set_nu(index_t i, double nu);
202
203
204			/**
205			* \brief Specifies a user vector where the centroids of the Laguerre
206			* cells will be stored after computing transport.
207			* \param[out] x a vector of nb points * dimension() doubles.
208			*/
209			void set_Laguerre_centroids(double* x) {
210		✗	Laguerre_centroids_ = x;
211			}
212
213			/**
214			* \brief Sets the maximum error.
215			* \param eps acceptable relative deviation for the measure of a
216			* Voronoi cell.
217			*/
218			void set_epsilon(double eps) {
219		✗	epsilon_ = eps;
220			}
221
222
223			/**
224			* \brief Sets the tolerance for linear solve.
225			* \param[in] eps the maximum value of
226			* \f$ \| Ax - b \| / \| b \| \f$b
227			*/
228			void set_linsolve_epsilon(double eps) {
229			linsolve_epsilon_ = eps;
230			}
231
232			/**
233			* \brief Sets the maximum number of iterations
234			* for linear solve.
235			* \param[in] maxiter the maximum number of iterations.
236			*/
237			void set_linsolve_maxiter(index_t maxiter) {
238			linsolve_maxiter_ = maxiter;
239			}
240
241			/**
242			* \brief Sets the maximum number of line search iterations.
243			* \param[in] maxiter the maximum number of step length reduction
244			* for line search.
245			*/
246			void set_linesearch_maxiter(index_t maxiter) {
247			linesearch_maxiter_ = maxiter;
248			}
249
250			/**
251			* \brief Sets the number of steplength reductions to be
252			* done at the first iteration.
253			* \param[in] init_iter the number of steplength reductions to
254			* be apllied at the first iteration. If left 0, start with
255			* Newton step at each iteration, else do tentative steplength
256			* prediction.
257			*/
258			void set_linesearch_init_iter(index_t init_iter) {
259			linesearch_init_iter_ = init_iter;
260			}
261
262			/**
263			* \brief Sets the weight of the regularization term.
264			* \details The regularization term (norm of the weight vector) cancels
265			* the translational degree of freedom of the weights.
266			* \param[in] eps_reg the weight of the regularization term. Use 0.0 for
267			* no regularization.
268			*/
269			void set_regularization(double eps_reg) {
270		✗	epsilon_regularization_ = eps_reg;
271			}
272
273
274			/**
275			* \brief Specifies whether a direct solver should be used.
276			* \param[in] solver one of OT_PRECG (default), OT_SUPERLU, OT_CHOLMOD.
277			* \details The direct solvers (OT_SUPERLU, OT_CHOLMOD) are recommended only for
278			* surfacic data, since the sparse factors become not so sparse when
279			* volumetric meshes are considered.
280			*/
281			void set_linear_solver(OTLinearSolver solver) {
282		✗	linear_solver_ = solver;
283		✗	}
284
285			/**
286			* \brief Computes the weights that realize the optimal
287			* transport map between the source mesh and the target
288			* pointset.
289			* \param[in] max_iterations maximum number of solver iterations.
290			*/
291			void optimize(index_t max_iterations);
292
293			/**
294			* \brief Enable/disable messages during optimization.
295			* \param[in] x true if messages should be displayed, false
296			* otherwise.
297			*/
298			void set_verbose(bool x) {
299		✗	verbose_ = x;
300			}
301
302			/**
303			* \brief Computes the weights that realize the optimal
304			* transport map between the source mesh and the target
305			* pointset.
306			* \details The algorithm is described in http://arxiv.org/abs/1603.05579
307			* Kitawaga, Merigot, Thibert, A Newton Algorithm for semi-discrete OT.
308			* \param[in] max_iterations maximum number of solver iterations.
309			* \param[in] n number of weights to optimize, used in hierarchical
310			* mode. If zero, optimizes all the weights.
311			*/
312			void optimize_full_Newton(index_t max_iterations, index_t n=0);
313
314			/**
315			* \brief Optimizes one level of the multilevel algorithm.
316			* \details The function supposes that the sequence [0,b)
317			* has been previously optimized. It is used to initialize
318			* the sequence [b,e). The whole sequence [0,e) is then
319			* optimized.
320			* \param[in] b index fo the first point in the level
321			* \param[in] e one position past the last index of the level
322			* \param[in] max_iterations maximum number of iterations
323			*/
324			void optimize_level(index_t b, index_t e, index_t max_iterations);
325
326			/**
327			* \brief Multi-level optimization.
328			* \details The points specified by set_points() need to have
329			* a hierarchical structure. They can be constructed by
330			* compute_hierarchical_sampling().
331			* \param[in] levels sample indices that correspond to level l are
332			* in the range levels[l] (included) ... levels[l+1] (excluded)
333			* \param[in] max_iterations maximum number of iterations
334			* \see compute_hierarchical_sampling()
335			*/
336			void optimize_levels(
337			const vector<index_t>& levels, index_t max_iterations
338			);
339
340			/**
341			* \brief Gets the number of points.
342			* \return The number of points, that was previously defined
343			* by set_points()
344			*/
345			index_t nb_points() const {
346			return weights_.size();
347			}
348
349			/**
350			* \brief Gets a point.
351			* \param[in] i index of the point
352			* \return a const pointer to the coordinates of the
353			* (dimension()+1)d point \p i
354			*/
355			const double* point_ptr(index_t i) const {
356			geo_debug_assert(i < (nb_points() + nb_air_particles()));
357		✗	return &(points_dimp1_[dimp1_ * i]);
358			}
359
360			/**
361			* \brief Gets weight of a point.
362			* \param[in] i index of the point
363			* \return the weight that was computed for point \p i
364			*/
365			double weight(index_t i) const {
366		✗	return weights_[i];
367			}
368
369			/**
370			* \brief Sets a weight of a point.
371			* \param[in] i index of the point
372			* \param[in] val new value of the weight
373			*/
374			void set_weight(index_t i, double val) {
375			weights_[i] = val;
376			}
377
378			/**
379			* \brief Gets the d+1-th coordinate of the embedding for a point.
380			* \param[in] i index of the point
381			* \return the d+1-th coordinate that was computed for point \p i
382			*/
383			double potential(index_t i) const {
384			return points_dimp1_[dimp1_*i + dimension_];
385			}
386
387			/**
388			* \brief Callback for the numerical solver.
389			* \details Evaluates the objective function and its gradient.
390			* \param[in] n number of variables
391			* \param[in] x current value of the variables
392			* \param[out] f current value of the objective function
393			* \param[out] g gradient of the objective function
394			*/
395			static void funcgrad_CB(
396			index_t n, double* x, double& f, double* g
397			);
398
399			/**
400			* \brief Callback for the numerical solver.
401			* \param[in] n number of variables
402			* \param[in] x current value of the variables
403			* \param[in] f current value of the objective function
404			* \param[in] g gradient of the objective function
405			* \param[in] gnorm norm of the gradient of the objective function
406			*/
407			static void newiteration_CB(
408			index_t n, const double* x, double f, const double* g, double gnorm
409			);
410
411			/**
412			* \brief Gets the restricted Voronoi diagram.
413			* \return a pointer to the restricted Voronoi diagram
414			*/
415			RestrictedVoronoiDiagram* RVD() {
416			return RVD_;
417			}
418
419			/**
420			* \brief Sets whether the restricted Voronoi diagram at
421			* each iteration should be saved.
422			* \details If flag is set, then each iteration is saved
423			* in file "RVD_nnn.geogram".
424			* \param[in] x true if each iteration should be saved,
425			* false otherwise.
426			* \param[in] show_RVD_seed if true, the seed associated
427			* with each restricted Voronoi cell is connected to it
428			* \param[in] last_iter_only if true, only the last iteration
429			* is saved
430			*/
431			void set_save_RVD_iter(
432			bool x,
433			bool show_RVD_seed = false,
434			bool last_iter_only = false
435			) {
436			if(last_iter_only) {
437			save_RVD_iter_ = false;
438			save_RVD_last_iter_ = true;
439			} else {
440			save_RVD_iter_ = x;
441			}
442			show_RVD_seed_ = show_RVD_seed;
443			}
444
445			/**
446			* \brief Computes a mesh with the restricted Voronoi diagram.
447			* \param[out] M a reference to the computed restricted Voronoi diagram.
448			*/
449			virtual void get_RVD(Mesh& M) = 0;
450
451			/**
452			* \brief Computes the centroids of the Laguerre cells.
453			* \param[out] centroids a pointer to the dimension()*nb_points
454			* coordinates of the centroids.
455			*/
456			virtual void compute_Laguerre_centroids(double* centroids) = 0;
457
458			/**
459			* \brief Updates the sparsity pattern of the Hessian right after
460			* a new Laguerre diagram was computed.
461			*/
462			void update_sparsity_pattern();
463
464			/**
465			* \brief Starts a new linear system.
466			* \param[in] n the dimension of the system
467			* \param[in] x pointer to a contiguous array of \p n doubles,
468			* where the solution will be stored.
469			*/
470			void new_linear_system(index_t n, double* x);
471
472			/**
473			* \brief Adds a coefficient to the matrix of the system.
474			* \param[in] i , j the indices of the coefficient
475			* \param[in] a the value to be added to the coefficient
476			*/
477		✗	void add_ij_coefficient(index_t i, index_t j, double a) {
478		✗	if(!user_H_g_) {
479		✗	nlAddIJCoefficient(i,j,a);
480			} else {
481		✗	if(user_H_ != nullptr) {
482			geo_debug_assert(user_H_->type == NL_MATRIX_SPARSE_DYNAMIC);
483		✗	nlSparseMatrixAdd((NLSparseMatrix*)user_H_, i, j, a);
484			}
485			}
486		✗	}
487
488			/**
489			* \brief Adds a coefficient to the right hand side.
490			* \param[in] i the index of the coefficient
491			* \param[in] a the value to be added to the coefficient
492			*/
493			void add_i_right_hand_side(index_t i, double a) {
494		✗	if(!user_H_g_) {
495		✗	nlAddIRightHandSide(i,a);
496			}
497			}
498
499			/**
500			* \brief Solves a linear system.
501			* \details The solution is stored in the vector that
502			* was previously specified to new_linear_system().
503			*/
504			void solve_linear_system();
505
506			/**
507			* \brief Sets the initial value of the weight associated
508			* with one of the points.
509			* \param[in] i index of the point, in 0..nb_points-1, where
510			* np_points corresponds to the parameter of set_points.
511			* \param[in] w the value of the weight.
512			*/
513			void set_initial_weight(index_t i, double w) {
514		✗	weights_[i] = w;
515			}
516
517			/**
518			* \brief Gets the total mass of the domain.
519			* \return the total mass.
520			*/
521			double total_mass() const {
522			return total_mass_;
523			}
524
525			/**
526			* \brief Computes the P1 Laplacian of the Laguerre cells.
527			* \param[in] Omega the domain, either a surfacic or a
528			* volumetric mesh.
529			* \param[in] weights the weights of the Laguerre diagram.
530			* \param[out] Laplacian P1 Laplacian of the Laguerre diagram or nullptr if
531			* not needed.
532			* \param[out] measures optional measures the measures of
533			* all Laguerre cells, or nullptr if not needed.
534			*/
535			void compute_P1_Laplacian(
536			const double* weights, NLMatrix Laplacian, double* measures
537			);
538
539			protected:
540
541			/**
542			* \brief Gets the mass of the Dirac associated with point p.
543			* \return the desired mass at point p.
544			*/
545			double nu(index_t p) const {
546		✗	return nu_.size() == 0 ? constant_nu_ : nu_[p];
547			}
548
549
550			/**
551			* \brief Callback for the numerical solver.
552			*/
553			virtual void newiteration();
554
555			/**
556			* \brief Saves the RVD at each iteration if
557			* specified on command line (just for debugging/
558			* explaining the algorithm).
559			* \param[in] id index to be used for the file, that
560			* will be named RVD_id.meshb
561			*/
562			void save_RVD(index_t id);
563
564			/**
565			* \brief Computes the objective function and its gradient.
566			* \param[in] n number of variables
567			* \param[in] w current value of the variables
568			* \param[out] f current value of the objective function
569			* \param[out] g gradient of the objective function
570			*/
571			void funcgrad(index_t n, double* w, double& f, double* g);
572
573			/**
574			* \brief Calls the callback for each intersection between a
575			* Laguerre cell and a simplex of the background mesh.
576			*/
577			virtual void call_callback_on_RVD() = 0;
578
579			/**
580			* \brief Computes the objective function, its gradient and its Hessian.
581			* \details Gradient and Hessian are used to solve a Newton
582			* step H p = -g
583			* \param[in] n number of variables
584			* \param[in] w current value of the variables
585			* \param[out] f current value of the objective function
586			* \param[out] g gradient of the objective function
587			*/
588			void eval_func_grad_Hessian(
589			index_t n, const double* w,
590			double& f, double* g
591			);
592
593			/**
594			* \brief Computes the stopping criterion of the solver.
595			* \details The stopping criterion is determined from
596			* the user-specified epsilon, number of samples and
597			* target measure of a cell (constant_nu_).
598			* \param n number of samples
599			* \return the gradient threshold
600			* \see set_epsilon()
601			*/
602			double gradient_threshold(index_t n) const {
603		✗	return ::sqrt(double(n) * geo_sqr(epsilon_ * constant_nu_));
604			}
605
606			public:
607
608			/**
609			* \brief Base class for the callbacks executed for each intersection
610			* between a Laguerre cell and a simplex of the background mesh.
611			*/
612			class Callback {
613			public:
614			/**
615			* \brief Callback constructor.
616			* \param[in] OTM a pointer to the OptimalTransportMap
617			*/
618		✗	Callback(
619			OptimalTransportMap* OTM
620		✗	) : OTM_(OTM),
621		✗	Newton_step_(false),
622		✗	eval_F_(false),
623		✗	n_(0),
624		✗	w_(nullptr),
625		✗	g_(nullptr),
626		✗	mg_(nullptr) {
627		✗	weighted_ =
628		✗	OTM->mesh().vertices.attributes().is_defined("weight");
629		✗	}
630
631			/**
632			* \brief Callback destructor.
633			*/
634			virtual ~Callback();
635
636			/**
637			* \brief Sets where centroids should be output.
638			* \details This computes mass times centroid. The mass can
639			* be retreived (and used to divide) from the gradient.
640			* \param[in] mg a pointer to the dimension()*nb_points coordinates
641			* of the centroids times the mass of the Laguerre cells
642			*/
643			void set_Laguerre_centroids(double* mg) {
644		✗	mg_ = mg;
645		✗	}
646
647			/**
648			* \brief Tests whether Laguerre centroids should be computed.
649			* \retval true if Laguerre centroids should be computed.
650			* \retval false otherwise.
651			*/
652			bool has_Laguerre_centroids() const {
653		✗	return (mg_ != nullptr);
654			}
655
656			/**
657			* \brief Gets a pointer to the Laguerre centroids.
658			* \return a pointer to nb vertices * dimension doubles
659			* with the centroids of the Laguerre cells times the mass
660			* of the Laguerre cells.
661			*/
662			double* Laguerre_centroids() {
663			return mg_;
664			}
665
666			/**
667			* \brief Sets the weight vector
668			* \param[in] w a const pointer to the weight vector.
669			* \param[in] n the number of weights in the weight vector.
670			*/
671			void set_w(const double* w, index_t n) {
672		✗	w_ = w;
673		✗	n_ = n;
674			}
675
676			/**
677			* \brief Specifies whether current step is a Newton step.
678			* \details If it is a Newton step, then the Hessian is
679			* computed.
680			* \param[in] Newton true if the current step is a Newton
681			* step, false otherwise.
682			*/
683			void set_Newton_step(bool Newton) {
684		✗	Newton_step_ = Newton;
685			}
686
687			/**
688			* \brief Tests whether the current step is a Newton step.
689			* \retval true if the current step is a Newton step.
690			* \retval false otherwise.
691			*/
692			bool is_Newton_step() const {
693		✗	return Newton_step_;
694			}
695
696			/**
697			* \brief Specifies the number of threads.
698			* \details This allocates one function value per thread.
699			* \param[in] nb the number of threads.
700			*/
701			void set_nb_threads(index_t nb) {
702		✗	funcval_.assign(nb, 0.0);
703			}
704
705			/**
706			* \brief Specifies where the gradient should be stored.
707			* \param[in] g a pointer to an array of nb points doubles.
708			*/
709			void set_g(double* g) {
710		✗	g_ = g;
711			}
712
713			/**
714			* \brief Specifies whether the objective function should
715			* be evaluated.
716			* \details The Newton solver does not need evaluating the
717			* objective function, only the BFGS solver needs it.
718			* \param[in] x true if the objective function should be
719			* evaluated, false otherwise. Default is false.
720			*/
721			void set_eval_F(bool x) {
722		✗	eval_F_ = x;
723			}
724
725			/**
726			* \brief Gets the computed value of the objective function.
727			* \details This sums the contributions of all threads.
728			* \retval the value of the objective function.
729			*/
730			double funcval() const {
731			double result = 0.0;
732		✗	FOR(i,funcval_.size()) {
733		✗	result += funcval_[i];
734			}
735			return result;
736			}
737
738			protected:
739			OptimalTransportMap* OTM_;
740			bool weighted_;
741			bool Newton_step_;
742			bool eval_F_;
743			vector<double> funcval_;
744			index_t n_;
745			const double* w_;
746			double* g_;
747			double* mg_;
748			};
749
750			protected:
751			static OptimalTransportMap* instance_;
752			index_t dimension_;
753			index_t dimp1_; /**< \brief dimension_ + 1 */
754			Mesh* mesh_;
755			Delaunay_var delaunay_;
756			RestrictedVoronoiDiagram_var RVD_;
757			vector<double> points_dimp1_;
758			vector<double> weights_;
759			double total_mass_;
760			double constant_nu_; /**< \brief Value of one of the Diracs if cte. */
761			vector<double> nu_; /**< \brief Value of all the Diracs. */
762			double epsilon_;
763			/**< \brief Acceptable relative deviation for the measure of a cell */
764			index_t current_call_iter_;
765
766			Callback* callback_;
767
768			std::string last_stats_;
769			bool pretty_log_;
770			index_t level_;
771
772			bool save_RVD_iter_;
773			bool save_RVD_last_iter_;
774			bool show_RVD_seed_;
775			index_t current_iter_;
776			bool newton_;
777			bool verbose_;
778
779			/**
780			* \brief Add a regularization term to remove
781			* translational degree of freedom for the
782			* weights.
783			*/
784			double epsilon_regularization_;
785
786			/**
787			* \brief Number of empty cells in last iteration.
788			*/
789			index_t nbZ_;
790
791			/**
792			* \brief Norm of the gradient in last iteration.
793			*/
794			double g_norm_;
795
796			/**
797			* \brief Measure of the smallest Laguerre cell.
798			*/
799			double measure_of_smallest_cell_;
800
801			/**
802			* \brief True if w did not change, thus there is
803			* no need to recompute the power diagram.
804			*/
805			bool w_did_not_change_;
806
807			/** \brief If user-specified, then Laguerre centroids are output here */
808			double* Laguerre_centroids_;
809
810			/** \brief maximum value of \f$ \| Ax - b \| / \| b \| \f$ */
811			double linsolve_epsilon_;
812
813			/** \brief maximum number of iterations for linear solve */
814			index_t linsolve_maxiter_;
815
816			/** \brief maximum number of steplength divisions */
817			index_t linesearch_maxiter_;
818
819			/** \brief starting number of steplength divisions */
820			index_t linesearch_init_iter_;
821
822			/** \brief one of OT_PRECG, OT_SUPERLU, OT_CHOLMOD. */
823			OTLinearSolver linear_solver_;
824
825			/** \brief if set, pointer to the air particles. */
826			const double* air_particles_;
827
828			/** \brief if non-zero, number of air particles. */
829			index_t nb_air_particles_;
830
831			/**
832			* \brief Number of doubles between two consecutive
833			* air particles in air_particles_.
834			*/
835			index_t air_particles_stride_;
836
837			/**
838			* \brief The fraction of the total mass occupied by air.
839			*/
840			double air_fraction_;
841
842			/**
843			* \brief Enabled if air fraction is specified without any
844			* air particles.
845			*/
846			bool clip_by_balls_;
847
848
849			/**
850			* \brief True if class is just used by user to compute
851			* Hessian and gradient instead of doing full computation.
852			*/
853			bool user_H_g_;
854
855			/**
856			* \brief User-defined Hessian matrix.
857			*/
858			NLMatrix user_H_;
859			};
860
861			}
862
863			#endif
864
865			#endif
866