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			#include <exploragram/optimal_transport/optimal_transport.h>
41			#include <exploragram/optimal_transport/linear_least_squares.h>
42
43			#include <geogram/mesh/mesh_io.h>
44			#include <geogram/mesh/mesh_reorder.h>
45			#include <geogram/mesh/mesh_geometry.h>
46			#include <geogram/mesh/mesh_AABB.h>
47
48			#include <geogram/voronoi/CVT.h>
49			#include <geogram/voronoi/generic_RVD_vertex.h>
50			#include <geogram/voronoi/RVD_callback.h>
51			#include <geogram/voronoi/generic_RVD_cell.h>
52			#include <geogram/delaunay/delaunay_3d.h>
53
54			#include <geogram/points/nn_search.h>
55			#include <geogram/numerics/optimizer.h>
56			#include <geogram/numerics/lbfgs_optimizers.h>
57
58			#include <geogram/basic/stopwatch.h>
59			#include <geogram/basic/file_system.h>
60			#include <geogram/basic/process.h>
61			#include <geogram/basic/progress.h>
62			#include <geogram/basic/permutation.h>
63			#include <geogram/basic/command_line.h>
64
65			#include <geogram/NL/nl.h>
66
67			#include <geogram/bibliography/bibliography.h>
68
69
70			namespace GEO {
71
72			OptimalTransportMap* OptimalTransportMap::instance_ = nullptr;
73
74		✗	OptimalTransportMap::Callback::~Callback() {
75		✗	}
76
77		✗	OptimalTransportMap::OptimalTransportMap(
78			index_t dimension,
79			Mesh* mesh, const std::string& delaunay, bool BRIO
80		✗	) : mesh_(mesh) {
81
82		✗	geo_cite("DBLP:conf/compgeom/AurenhammerHA92");
83		✗	geo_cite("DBLP:journals/cgf/Merigot11");
84		✗	geo_cite("journals/M2AN/LevyNAL15");
85
86		✗	dimension_ = dimension;
87		✗	dimp1_ = dimension_+1;
88
89		✗	epsilon_regularization_ = 0.0;
90
91			// Mesh is supposed to be embedded in d+1 dim (with
92			// (d+1-th dimension set to zero).
93		✗	geo_assert(mesh->vertices.dimension() == dimp1_);
94
95			// Note: we represent power diagrams as d+1 Voronoi diagrams
96		✗	delaunay_ = Delaunay::create(coord_index_t(dimp1_), delaunay);
97
98		✗	RVD_ = RestrictedVoronoiDiagram::create(delaunay_, mesh_);
99		✗	RVD_->set_volumetric(true);
100		✗	RVD_->set_check_SR(true);
101		✗	RVD_->create_threads();
102
103			// No need to reorder vertices if BRIO is activated since
104			// vertices are then already reordered.
105		✗	if(BRIO) {
106		✗	RVD_->delaunay()->set_reorder(false);
107			}
108
109		✗	newton_ = false;
110		✗	verbose_ = true;
111
112		✗	instance_ = nullptr;
113		✗	constant_nu_ = 0.0;
114		✗	total_mass_ = 0.0;
115		✗	current_call_iter_ = 0;
116		✗	epsilon_ = 0.01;
117		✗	level_ = 0;
118
119		✗	save_RVD_iter_ = false;
120		✗	save_RVD_last_iter_ = false;
121		✗	show_RVD_seed_ = false;
122		✗	current_iter_ = 0;
123
124		✗	pretty_log_ = false; // CmdLine::get_arg_bool("log:pretty");
125
126		✗	w_did_not_change_ = false;
127		✗	measure_of_smallest_cell_ = 0.0;
128		✗	callback_ = nullptr;
129		✗	Laguerre_centroids_ = nullptr;
130
131		✗	linsolve_epsilon_ = 0.001;
132		✗	linsolve_maxiter_ = 1000;
133
134		✗	linesearch_maxiter_ = 100;
135		✗	linesearch_init_iter_ = 0;
136
137		✗	linear_solver_ = OT_PRECG;
138
139		✗	nb_air_particles_ = 0;
140		✗	air_particles_ = nullptr;
141		✗	air_particles_stride_ = 0;
142		✗	air_fraction_ = 0.0;
143
144		✗	clip_by_balls_ = false;
145
146		✗	user_H_g_ = false;
147		✗	user_H_ = nullptr;
148		✗	}
149
150		✗	OptimalTransportMap::~OptimalTransportMap() {
151		✗	delete callback_;
152		✗	callback_ = nullptr;
153		✗	}
154
155		✗	void OptimalTransportMap::set_points(
156			index_t nb_points, const double* points, index_t stride
157			) {
158
159		✗	if(stride == 0) {
160		✗	stride = dimension_;
161			}
162
163		✗	index_t nb_total = nb_points + nb_air_particles_;
164
165			// Note: we represent power diagrams as (d+1)-dim Voronoi diagrams.
166			// The target points are lifted to (d+1)-dim.
167		✗	points_dimp1_.resize(nb_total * dimp1_);
168		✗	for(index_t i = 0; i < nb_points; ++i) {
169		✗	double* p = &points_dimp1_[i*dimp1_];
170		✗	for(index_t c=0; c<dimension_; ++c) {
171		✗	p[c] = points[i*stride+c];
172			}
173		✗	p[dimension()] = 0.0; // Yes, dimension() and not dimension()-1
174			// (for instance, in 2d, x->0, y->1, W->2)
175			}
176		✗	for(index_t i=0; i<nb_air_particles_; ++i) {
177		✗	double* p = &points_dimp1_[(i + nb_points)*dimp1_];
178		✗	for(index_t c=0; c<dimension_; ++c) {
179		✗	p[c] = air_particles_[i*air_particles_stride_+c];
180			}
181		✗	p[dimension()] = 0.0; // Yes, dimension() and not dimension()-1
182			// (for instance, in 2d, x->0, y->1, W->2)
183			}
184		✗	weights_.assign(nb_points, 0);
185		✗	constant_nu_ = (1.0 - air_fraction_) * total_mass_ / double(nb_points);
186		✗	if(air_fraction_ != 0.0 && nb_air_particles_ == 0) {
187			// constant_nu_ = pi*R^2 -> R^2 = constant_nu_ / pi
188			// TODO: volumetric case.
189		✗	weights_.assign(nb_points, constant_nu_ / M_PI);
190			}
191		✗	}
192
193		✗	void OptimalTransportMap::set_nu(index_t i, double nu) {
194		✗	geo_debug_assert(i < weights_.size());
195		✗	if(nu_.size() != weights_.size()) {
196		✗	nu_.assign(weights_.size(), 0.0);
197			}
198		✗	nu_[i] = nu;
199		✗	}
200
201		✗	void OptimalTransportMap::optimize_full_Newton(
202			index_t max_iterations, index_t n
203			) {
204		✗	if(n == 0) {
205		✗	n = index_t(points_dimp1_.size() / dimp1_) - nb_air_particles_;
206			}
207
208		✗	vector<double> pk(n);
209		✗	vector<double> xk(n);
210		✗	vector<double> gk(n);
211		✗	double fk=0.0;
212		✗	bool converged = false;
213
214		✗	double epsilon0 = 0.0;
215		✗	w_did_not_change_ = false;
216
217		✗	if(max_iterations == 0) {
218		✗	if(Laguerre_centroids_ != nullptr) {
219		✗	callback_->set_Laguerre_centroids(Laguerre_centroids_);
220			}
221		✗	funcgrad(n,weights_.data(),fk,gk.data());
222		✗	if(Laguerre_centroids_ != nullptr) {
223		✗	callback_->set_Laguerre_centroids(nullptr);
224			}
225			}
226
227			// Inner iteration control for linesearch
228		✗	index_t first_inner_iter = linesearch_init_iter_;
229		✗	index_t inner_iter = first_inner_iter;
230		✗	bool use_inner_iter_prediction = (first_inner_iter != 0);
231
232		✗	for(index_t k=0; k<max_iterations; ++k) {
233		✗	if(verbose_) {
234		✗	std::cerr << "======= k = " << k << std::endl;
235			}
236		✗	xk=weights_;
237
238		✗	new_linear_system(n,pk.data());
239		✗	eval_func_grad_Hessian(n,xk.data(),fk,gk.data());
240
241		✗	if(k == 0) {
242		✗	newiteration();
243			}
244
245		✗	if(epsilon0 == 0.0) {
246		✗	epsilon0 = measure_of_smallest_cell_;
247		✗	FOR(i,n) {
248		✗	epsilon0 = std::min(epsilon0, nu(i));
249			}
250		✗	if(nb_air_particles_ != 0.0) {
251		✗	epsilon0 = std::min(epsilon0, air_fraction_ * total_mass_);
252			}
253		✗	epsilon0 = 0.5 * epsilon0;
254			}
255
256		✗	if(verbose_) {
257		✗	Logger::out("OTM") << " Solving linear system" << std::endl;
258			}
259
260		✗	if(nbZ_ != 0) {
261		✗	std::cerr << "There were empty cells !!!!!!" << std::endl;
262		✗	std::cerr << "FATAL error, exiting Newton" << std::endl;
263		✗	return;
264			}
265		✗	solve_linear_system();
266
267		✗	if(verbose_) {
268		✗	std::cerr << "Line search ..." << std::endl;
269			}
270
271		✗	w_did_not_change_ = false;
272
273		✗	double alphak = 1.0;
274		✗	double gknorm = g_norm_;
275
276		✗	if(first_inner_iter != 0) {
277		✗	alphak /= pow(2.0, double(first_inner_iter));
278			}
279
280		✗	for(inner_iter=first_inner_iter;
281		✗	inner_iter < linesearch_maxiter_; ++inner_iter
282			) {
283		✗	if(verbose_) {
284		✗	std::cerr << " inner iter = "
285		✗	<< inner_iter << std::endl;
286			}
287
288			// weights = xk + alphak pk
289		✗	for(index_t i=0; i<n; ++i) {
290		✗	weights_[i] = xk[i] + alphak * pk[i];
291			}
292
293			// Compute cell measures and nbZ.
294		✗	if(Laguerre_centroids_ != nullptr) {
295		✗	callback_->set_Laguerre_centroids(Laguerre_centroids_);
296			}
297
298		✗	funcgrad(n,weights_.data(),fk,gk.data());
299
300		✗	if(Laguerre_centroids_ != nullptr) {
301		✗	callback_->set_Laguerre_centroids(nullptr);
302			}
303
304		✗	if(verbose_) {
305		✗	std::cerr << "cell measure :"
306		✗	<< measure_of_smallest_cell_
307		✗	<< "(>=?)" << epsilon0 << std::endl;
308		✗	std::cerr << "gradient norm:"
309		✗	<< g_norm_ << "(<=?)"
310		✗	<< (1.0 - 0.5*alphak) * gknorm << std::endl;
311			}
312
313			// Condition to exit linesearch loop
314		✗	if(
315		✗	(measure_of_smallest_cell_ >= epsilon0) && (
316		✗	(air_fraction_ != 0.0) ||
317		✗	(g_norm_ <= (1.0 - 0.5*alphak) * gknorm)
318			)
319			) {
320			// Condition for global convergence
321		✗	if(g_norm_ < gradient_threshold(n)) {
322		✗	converged = true;
323			}
324		✗	break;
325			}
326			// Else we halve the step.
327		✗	alphak /= 2.0;
328			}
329
330		✗	if(use_inner_iter_prediction) {
331		✗	if(inner_iter <= 2) {
332		✗	first_inner_iter = 0;
333			} else {
334		✗	first_inner_iter = inner_iter / 2;
335			}
336			} else {
337		✗	first_inner_iter = 0;
338			}
339
340		✗	newiteration();
341		✗	if(converged) {
342		✗	break;
343			}
344			// No need to update the power diagram at next iteration,
345			// since we will evaluate the Hessian for the same weight
346			// vector.
347		✗	w_did_not_change_ = true;
348			}
349		✗	if(save_RVD_last_iter_) {
350		✗	save_RVD(current_iter_);
351			}
352		✗	}
353
354		✗	void OptimalTransportMap::optimize(index_t max_iterations) {
355
356		✗	index_t n = index_t(points_dimp1_.size() / dimp1_) - nb_air_particles_;
357
358			// Sanity check
359		✗	if(nu_.size() != 0) {
360		✗	double total_nu = 0.0;
361		✗	FOR(i,n) {
362		✗	total_nu += nu(i);
363			}
364		✗	if(verbose_) {
365		✗	std::cerr << "total nu=" << total_nu << std::endl;
366		✗	std::cerr << "total mass=" << total_mass_ << std::endl;
367			}
368		✗	if(::fabs(total_nu - total_mass_)/total_mass_ > 0.01) {
369		✗	Logger::warn("OTM")
370		✗	<< "Specified nu do not sum to domain measure"
371		✗	<< std::endl;
372		✗	Logger::warn("OTM")
373		✗	<< "rescaling..."
374		✗	<< std::endl;
375			}
376		✗	FOR(i,n) {
377		✗	set_nu(i, nu(i) * total_mass_ / total_nu);
378			}
379
380		✗	total_nu = 0.0;
381		✗	FOR(i,n) {
382		✗	total_nu += nu(i);
383			}
384		✗	if(::fabs(total_nu - total_mass_)/total_mass_ > 0.01) {
385		✗	Logger::warn("OTM")
386		✗	<< "Specified nu do not sum to domain measure"
387		✗	<< std::endl;
388		✗	return;
389			}
390			}
391
392		✗	if(newton_) {
393		✗	optimize_full_Newton(max_iterations);
394		✗	return;
395			}
396
397		✗	level_ = 0;
398		✗	index_t m = 7;
399		✗	Optimizer_var optimizer = Optimizer::create("HLBFGS");
400
401		✗	optimizer->set_epsg(gradient_threshold(n));
402		✗	optimizer->set_epsf(0.0);
403		✗	optimizer->set_epsx(0.0);
404
405		✗	optimizer->set_newiteration_callback(newiteration_CB);
406		✗	optimizer->set_funcgrad_callback(funcgrad_CB);
407
408		✗	optimizer->set_N(n);
409		✗	optimizer->set_M(m);
410		✗	optimizer->set_max_iter(max_iterations);
411		✗	instance_ = this;
412		✗	current_call_iter_ = 0;
413		✗	current_iter_ = 0;
414
415		✗	callback_->set_eval_F(true);
416		✗	optimizer->optimize(weights_.data());
417		✗	callback_->set_eval_F(false);
418
419		✗	instance_ = nullptr;
420			// To make sure everything is reset properly
421		✗	double dummy = 0;
422		✗	funcgrad(n, weights_.data(), dummy, nullptr);
423		✗	if(verbose_) {
424		✗	Logger::out("OTM")
425		✗	<< "Used " << current_call_iter_ << " iterations" << std::endl;
426			}
427		✗	if(save_RVD_last_iter_) {
428		✗	save_RVD(current_iter_);
429			}
430		✗	}
431
432
433		✗	void OptimalTransportMap::optimize_level(
434			index_t b, index_t e, index_t max_iterations
435			) {
436
437			// If this is not the first level, propagate the weights from
438			// the lower levels.
439		✗	if(b != 0) {
440
441			// Create a nearest neighbor search data structure
442			// and insert the [0..b) samples into it (they were
443			// initialized at previous calls).
444			NearestNeighborSearch_var NN =
445		✗	NearestNeighborSearch::create(coord_index_t(dimension()));
446
447		✗	NN->set_points(b, points_dimp1_.data(), dimp1_);
448		✗	index_t degree = 2; // CmdLine::get_arg_uint("fitting_degree");
449
450			// If degree \notin {1,2}, use weight of nearest sample
451		✗	if(degree < 1 || degree > 2) {
452		✗	for(index_t i = b; i < e; ++i) {
453		✗	weights_[i] =
454			weights_[
455		✗	NN->get_nearest_neighbor(&points_dimp1_[dimp1_ * i])
456		✗	];
457			}
458		✗	} else {
459
460			// If degree \in {1,2} use linear least squares to
461			// compute an estimate of the weight function.
462		✗	LinearLeastSquares LLS(degree);
463		✗	for(index_t i = b; i < e; ++i) {
464		✗	const index_t nb = 10 * degree;
465			index_t neighbor[100];
466			double dist[100];
467		✗	NN->get_nearest_neighbors(
468		✗	nb, &points_dimp1_[dimp1_ * i], neighbor, dist
469			);
470		✗	LLS.begin();
471		✗	for(index_t jj = 0; jj < nb; ++jj) {
472		✗	if(dist[jj] != 0.0) {
473		✗	index_t j = neighbor[jj];
474		✗	LLS.add_point(
475		✗	&points_dimp1_[dimp1_ * j], weights_[j]
476			);
477			}
478			}
479		✗	LLS.end();
480		✗	weights_[i] = LLS.eval(&points_dimp1_[dimp1_ * i]);
481			}
482			}
483		✗	}
484
485			// Optimize the weights associated with the sequence [0,e)
486		✗	index_t n = e;
487
488			// Important! constant_nu_ (target measure of a cell) needs
489			// to be updated, since it depends on the number of samples
490			// (that varies at each level).
491		✗	constant_nu_ = total_mass_ / double(n);
492
493		✗	if(newton_) {
494		✗	optimize_full_Newton(max_iterations, n);
495		✗	return;
496			}
497
498		✗	index_t m = 7;
499		✗	Optimizer_var optimizer = Optimizer::create("HLBFGS");
500
501		✗	optimizer->set_epsg(gradient_threshold(n));
502		✗	optimizer->set_epsf(0.0);
503		✗	optimizer->set_epsx(0.0);
504		✗	optimizer->set_newiteration_callback(newiteration_CB);
505		✗	optimizer->set_funcgrad_callback(funcgrad_CB);
506
507		✗	optimizer->set_N(n);
508		✗	optimizer->set_M(m);
509		✗	optimizer->set_max_iter(max_iterations);
510		✗	instance_ = this;
511		✗	current_call_iter_ = 0;
512		✗	callback_->set_eval_F(true);
513		✗	optimizer->optimize(weights_.data());
514		✗	callback_->set_eval_F(false);
515		✗	instance_ = nullptr;
516
517			// To make sure everything is reset properly
518		✗	double dummy = 0;
519		✗	funcgrad(n, weights_.data(), dummy, nullptr);
520		✗	}
521
522		✗	void OptimalTransportMap::optimize_levels(
523			const vector<index_t>& levels, index_t max_iterations
524			) {
525		✗	if(verbose_) {
526		✗	if(levels.size() > 2) {
527		✗	Logger::out("OTM") << "Using " << levels.size()-1
528		✗	<< " levels" << std::endl;
529			} else {
530		✗	Logger::out("OTM") << "Using 1 level" << std::endl;
531			}
532			}
533		✗	for(index_t l = 0; l + 1 < levels.size(); ++l) {
534		✗	level_ = l+1;
535		✗	index_t b = levels[l];
536		✗	index_t e = levels[l + 1];
537		✗	vector<index_t> brio_levels;
538		✗	for(index_t i=0; i<=l+1; ++i) {
539		✗	brio_levels.push_back(levels[i]);
540			}
541		✗	RVD_->delaunay()->set_BRIO_levels(brio_levels);
542		✗	optimize_level(b, e, max_iterations);
543		✗	}
544		✗	if(save_RVD_last_iter_) {
545		✗	save_RVD(current_iter_);
546			}
547		✗	}
548
549		✗	void OptimalTransportMap::funcgrad_CB(
550			index_t n, double* x, double& f, double* g
551			) {
552		✗	instance_->funcgrad(n, x, f, g);
553		✗	}
554
555		✗	void OptimalTransportMap::newiteration_CB(
556			index_t n, const double* x, double f, const double* g, double gnorm
557			) {
558		✗	geo_argused(n);
559		✗	geo_argused(x);
560		✗	geo_argused(f);
561		✗	geo_argused(g);
562		✗	geo_argused(gnorm);
563		✗	instance_->newiteration();
564		✗	}
565
566		✗	void OptimalTransportMap::newiteration() {
567			//xxx std::cerr << "newiteration" << std::endl;
568		✗	if(save_RVD_iter_) {
569		✗	if(verbose_) {
570		✗	std::cerr << " save iter" << std::endl;
571			}
572		✗	save_RVD(current_iter_);
573			}
574		✗	++current_iter_;
575		✗	}
576
577		✗	void OptimalTransportMap::save_RVD(index_t id) {
578		✗	Mesh RVD_mesh;
579		✗	get_RVD(RVD_mesh);
580		✗	MeshIOFlags flags;
581		✗	flags.set_attribute(MESH_CELL_REGION);
582		✗	flags.set_attribute(MESH_FACET_REGION);
583		✗	mesh_save(
584			RVD_mesh,
585		✗	"RVD_" + String::to_string(id) + ".geogram",
586			flags
587			);
588		✗	}
589
590		✗	void OptimalTransportMap::funcgrad(
591			index_t n, double* w, double& f, double* g
592			) {
593
594		✗	bool is_Newton_step = callback_->is_Newton_step();
595
596			// For now, always compute function and gradient
597		✗	bool update_fg = true;
598
599			// Delaunay triangulation is only updated if function and
600			// gradient is evaluated. If only Hessian needs to be evaluated,
601			// then it is at the same point as the latest function and gradient
602			// evaluation (see Yang Liu's CVT-Newton code).
603		✗	if(update_fg && !w_did_not_change_) {
604			// Step 1: determine the (dim+1)d embedding from the weights
605		✗	double W = 0.0;
606		✗	for(index_t p = 0; p < n; ++p) {
607		✗	W = std::max(W, w[p]);
608			}
609		✗	for(index_t p = 0; p < n; ++p) {
610			// Yes, dimension_ and not dimension_ -1,
611			// for instance in 2d, x->0, y->1, W->2
612		✗	points_dimp1_[dimp1_ * p + dimension_] = ::sqrt(W - w[p]);
613			}
614		✗	if(nb_air_particles_ != 0) {
615		✗	for(index_t p = 0; p < nb_air_particles_; ++p) {
616		✗	points_dimp1_[dimp1_ * (n + p) + dimension_] =
617		✗	::sqrt(W - 0.0);
618			}
619			}
620
621			// Step 2: compute function and gradient
622			{
623		✗	Stopwatch* SW = nullptr;
624		✗	if(newton_) {
625		✗	if(verbose_) {
626		✗	SW = new Stopwatch("Power diagram");
627		✗	Logger::out("OTM") << "In power diagram..."
628		✗	<< std::endl;
629			}
630			}
631		✗	delaunay_->set_vertices(
632		✗	(n + nb_air_particles_), points_dimp1_.data()
633			);
634		✗	if(verbose_ && newton_) {
635		✗	delete SW;
636			}
637			}
638			}
639
640		✗	if(is_Newton_step) {
641		✗	update_sparsity_pattern();
642			}
643
644		✗	if(g == nullptr) {
645		✗	if(pretty_log_) {
646		✗	CmdLine::ui_clear_line();
647		✗	CmdLine::ui_message(last_stats_ + "\n");
648			}
649		✗	return;
650			}
651
652		✗	if(update_fg) {
653		✗	f = 0.0;
654		✗	for(index_t p = 0; p < n; ++p) {
655		✗	g[p] = 0.0;
656			}
657			}
658
659		✗	callback_->set_w(w,n);
660		✗	callback_->set_g(g);
661		✗	callback_->set_nb_threads(Process::maximum_concurrent_threads());
662
663		✗	if(callback_->has_Laguerre_centroids()) {
664		✗	Memory::clear(
665		✗	callback_->Laguerre_centroids(),
666		✗	nb_points()*sizeof(double)*dimension()
667			);
668			}
669
670			{
671		✗	Stopwatch* W = nullptr;
672		✗	if(verbose_ && newton_) {
673		✗	W = new Stopwatch("RVD");
674		✗	Logger::out("OTM") << "In RVD (funcgrad)..." << std::endl;
675			}
676		✗	call_callback_on_RVD();
677		✗	if(verbose_ && newton_) {
678		✗	delete W;
679			}
680			}
681		✗	f = callback_->funcval();
682
683		✗	if(callback_->has_Laguerre_centroids()) {
684		✗	for(index_t v=0; v<nb_points(); ++v) {
685		✗	for(index_t c=0; c<dimension_; ++c) {
686		✗	callback_->Laguerre_centroids()[dimension_*v+c] /= g[v];
687			}
688			}
689			}
690
691		✗	index_t nb_empty_cells = 0;
692		✗	if(update_fg) {
693		✗	measure_of_smallest_cell_ = Numeric::max_float64();
694		✗	for(index_t i=0; i<n; ++i) {
695		✗	measure_of_smallest_cell_ =
696		✗	std::min(measure_of_smallest_cell_, g[i]);
697		✗	if(g[i] == 0.0) {
698		✗	++nb_empty_cells;
699			}
700			}
701		✗	if(air_fraction_ != 0.0) {
702		✗	double air_mass = total_mass_;
703		✗	for(index_t i=0; i<n; ++i) {
704		✗	air_mass -= g[i];
705			}
706		✗	if(fabs(air_mass) < 1e-30) {
707		✗	++nb_empty_cells;
708			}
709		✗	measure_of_smallest_cell_ =
710		✗	std::min(measure_of_smallest_cell_, air_mass);
711			}
712			}
713
714		✗	if(update_fg) {
715		✗	for(index_t p = 0; p < n; ++p) {
716		✗	f += nu(p) * w[p];
717			// Note: we minimize -f instead of maximizing f,
718			// therefore, in the paper:
719			// g[p] = lambda_p - mesure(power cell associated with p)
720			//
721			// What is programmed:
722			// g[p] = mesure(power cell associated with p) - lambda_p
723		✗	g[p] -= nu(p);
724
725		✗	if(is_Newton_step) {
726			// Newton step: solve H deltax = -g
727			// (note the minus sign on the right hand side)
728			// g[p] -= lamnda_p_ -> RHS[p] += constant_nu_
729		✗	add_i_right_hand_side(p, nu(p));
730			}
731			}
732			}
733
734		✗	double max_diff = 0.0;
735		✗	double avg_diff = 0.0;
736
737		✗	for(index_t p = 0; p < n; ++p) {
738		✗	double cur_diff = ::fabs(g[p]);
739		✗	max_diff = std::max(max_diff, cur_diff);
740		✗	avg_diff += cur_diff / double(n);
741			}
742
743
744			// Regularisation: minimize the squared norm of the weight
745			// vector to remove a translational degree of freedom.
746			// It seems to make the overall convergence slower
747			// (but this may be due to a wrong
748			// scaling between the different levels, to be investigated...)
749		✗	if(epsilon_regularization_ != 0.0) {
750		✗	if(update_fg) {
751		✗	for(index_t p = 0; p < n; ++p) {
752		✗	f += 0.5 * epsilon_regularization_ * nu(p) * w[p]*w[p];
753		✗	g[p] += epsilon_regularization_ * nu(p) * w[p];
754			}
755			}
756		✗	if(is_Newton_step) {
757		✗	for(index_t p = 0; p < n; ++p) {
758		✗	add_ij_coefficient(
759		✗	p,p,epsilon_regularization_*nu(p)
760			);
761		✗	add_i_right_hand_side(
762		✗	p,-epsilon_regularization_*nu(p)*w[p]
763			);
764			}
765			}
766			}
767
768
769		✗	double gNorm = 0.0;
770		✗	for(index_t i = 0; i < n; ++i) {
771		✗	gNorm += geo_sqr(g[i]);
772			}
773		✗	gNorm = ::sqrt(gNorm);
774
775		✗	nbZ_ = nb_empty_cells;
776
777		✗	std::ostringstream str;
778		✗	if(pretty_log_) {
779		✗	if(level_ == 0) {
780		✗	str << "o-[OTM ] " ;
781			} else {
782		✗	str << "o-[OTM Lvl." << level_ << " ] " ;
783			}
784			} else {
785		✗	if(level_ == 0) {
786		✗	str << " OTM : " ;
787			} else {
788		✗	str << " OTM Lvl." << level_ << ": " ;
789			}
790			}
791
792		✗	double scl = 100.0 / constant_nu_;
793
794		✗	str << "iter=" << current_call_iter_
795		✗	<< " nbZ=" << nb_empty_cells
796			// << " f=" << f
797		✗	<< " avg_diff=" << (avg_diff * scl) << "%"
798		✗	<< " max_diff=" << (max_diff * scl) << "%"
799		✗	<< " g=" << gNorm
800			// << " f=" << f
801		✗	<< " threshold=" << gradient_threshold(n);
802		✗	last_stats_ = str.str();
803
804		✗	g_norm_ = gNorm;
805
806			// "custom task progress" (clears the standard message
807			// and replaces it with another one).
808		✗	if(verbose_) {
809		✗	if(pretty_log_) {
810		✗	if(current_call_iter_ != 0) {
811		✗	CmdLine::ui_clear_line();
812			}
813		✗	CmdLine::ui_message(str.str());
814			} else {
815			// str << " f=" << f;
816		✗	CmdLine::ui_message(str.str() + "\n");
817			}
818			}
819		✗	++current_call_iter_;
820		✗	}
821
822		✗	void OptimalTransportMap::eval_func_grad_Hessian(
823			index_t n, const double* w, double& f, double* g
824			) {
825		✗	callback_->set_Newton_step(true);
826		✗	funcgrad(n,const_cast<double*>(w),f,g);
827		✗	callback_->set_Newton_step(false);
828		✗	}
829
830			/************************************************************/
831
832			// TODO: in the Euler code, see if we do
833			// not have duplicated computations, i.e.
834			// - Power diagrams when leaving and entering iteration ?
835			// - Centroids: do we restart a RVD computation ?
836			// TODO: are we obliged to create/destroy OpenNL context for each system ?
837			// TODO: we could have an OpenNL buffer for the RHS of the Newton solve ?
838			// Not really, because this is *minus* the gradient.
839
840		✗	void OptimalTransportMap::update_sparsity_pattern() {
841			// Does nothing for now,
842			// (we let OpenNL discover the sparsity pattern)
843			// Tryed smarter things, but was not faster...
844		✗	}
845
846		✗	void OptimalTransportMap::new_linear_system(index_t n, double* x) {
847		✗	nlNewContext();
848
849		✗	if(linear_solver_ != OT_PRECG) {
850		✗	if(linear_solver_ == OT_SUPERLU) {
851		✗	if(nlInitExtension("SUPERLU") != NL_TRUE) {
852		✗	linear_solver_ = OT_PRECG;
853			}
854		✗	} else if(linear_solver_ == OT_CHOLMOD) {
855		✗	if(nlInitExtension("CHOLMOD") != NL_TRUE) {
856		✗	linear_solver_ = OT_PRECG;
857			}
858			}
859		✗	if(linear_solver_ == OT_PRECG) {
860		✗	Logger::warn("OTM") << "Could not initialize OpenNL extension"
861		✗	<< std::endl;
862		✗	Logger::warn("OTM") << "Falling back to conjugate gradient"
863		✗	<< std::endl;
864			}
865			}
866
867		✗	if(verbose_) {
868		✗	nlEnable(NL_VERBOSE);
869			}
870
871		✗	nlSolverParameteri(NL_NB_VARIABLES, NLint(n));
872		✗	switch(linear_solver_) {
873		✗	case OT_PRECG:
874		✗	nlSolverParameteri(NL_SOLVER, NL_CG);
875		✗	nlSolverParameteri(NL_PRECONDITIONER, NL_PRECOND_JACOBI);
876		✗	nlSolverParameteri(NL_SYMMETRIC, NL_TRUE);
877		✗	nlSolverParameterd(NL_THRESHOLD, linsolve_epsilon_);
878		✗	nlSolverParameteri(NL_MAX_ITERATIONS, NLint(linsolve_maxiter_));
879
880		✗	break;
881		✗	case OT_SUPERLU:
882		✗	nlSolverParameteri(NL_SOLVER, NL_PERM_SUPERLU_EXT);
883		✗	break;
884		✗	case OT_CHOLMOD:
885		✗	nlSolverParameteri(NL_SOLVER, NL_CHOLMOD_EXT);
886		✗	break;
887			}
888		✗	nlEnable(NL_VARIABLES_BUFFER);
889		✗	nlEnable(NL_NO_VARIABLES_INDIRECTION);
890		✗	nlBegin(NL_SYSTEM);
891		✗	nlBindBuffer(NL_VARIABLES_BUFFER, 0, x, NLuint(sizeof(double)));
892		✗	nlBegin(NL_MATRIX);
893		✗	}
894
895		✗	void OptimalTransportMap::solve_linear_system() {
896		✗	nlEnd(NL_MATRIX);
897		✗	nlEnd(NL_SYSTEM);
898		✗	nlSolve();
899		✗	if(verbose_) {
900			int used_iters;
901			double elapsed_time;
902			double gflops;
903			double error;
904		✗	nlGetIntegerv(NL_USED_ITERATIONS, &used_iters);
905		✗	nlGetDoublev(NL_ELAPSED_TIME, &elapsed_time);
906		✗	nlGetDoublev(NL_GFLOPS, &gflops);
907		✗	nlGetDoublev(NL_ERROR, &error);
908		✗	std::cerr << " "
909		✗	<< used_iters << " iters in "
910		✗	<< elapsed_time << " seconds "
911		✗	<< gflops << " GFlop/s"
912		✗	<< " ||Ax-b||/||b||="
913		✗	<< error
914		✗	<< std::endl;
915			}
916		✗	nlDeleteContext(nlGetCurrent());
917		✗	}
918
919		✗	void OptimalTransportMap::compute_P1_Laplacian(
920			const double* w, NLMatrix Laplacian, double* measures
921			) {
922		✗	index_t n = index_t(points_dimp1_.size() / dimp1_) - nb_air_particles_;
923
924		✗	if(measures != nullptr) {
925		✗	Memory::clear(measures, n*sizeof(double));
926			}
927
928		✗	user_H_g_ = true;
929		✗	user_H_ = Laplacian;
930		✗	callback_->set_Newton_step(true);
931
932			// Step 1: determine the (dim+1)d embedding from the weights
933		✗	double W = 0.0;
934		✗	for(index_t p = 0; p < n; ++p) {
935		✗	W = std::max(W, w[p]);
936			}
937		✗	for(index_t p = 0; p < n; ++p) {
938			// Yes, dimension_ and not dimension_ -1,
939			// for instance in 2d, x->0, y->1, W->2
940		✗	points_dimp1_[dimp1_ * p + dimension_] = ::sqrt(W - w[p]);
941			}
942
943		✗	if(nb_air_particles_ != 0) {
944		✗	for(index_t p = 0; p < nb_air_particles_; ++p) {
945		✗	points_dimp1_[dimp1_ * (n + p) + dimension_] = ::sqrt(W - 0.0);
946			}
947			}
948
949			// Step 2: compute Laplacian and cell measures.
950		✗	delaunay_->set_vertices((n + nb_air_particles_), points_dimp1_.data());
951		✗	callback_->set_w(w,n);
952		✗	callback_->set_g(measures);
953		✗	callback_->set_nb_threads(Process::maximum_concurrent_threads());
954		✗	call_callback_on_RVD();
955
956		✗	callback_->set_Newton_step(false);
957		✗	user_H_g_ = false;
958		✗	user_H_ = nullptr;
959		✗	}
960
961
962
963			/**********************************************************************/
964			}
965