GCC Code Coverage Report

Line	Branch	Exec	Source
1
2			#ifdef GEOGRAM_WITH_HLBFGS
3
4			#include <geogram/basic/command_line.h>
5			#include <geogram/basic/command_line_args.h>
6			#include <geogram/numerics/optimizer.h>
7			#include <geogram/numerics/lbfgs_optimizers.h>
8			#include <geogram/third_party/HLBFGS/Lite_Sparse_Matrix.h>
9			#include <geogram/third_party/HLBFGS/HLBFGS.h>
10			#include <iostream>
11
12			namespace {
13			Lite_Sparse_Matrix* m_sparse_matrix = nullptr;
14
15			/********************** Initial HLBFGS API ***********************************/
16
17		✗	void evalfunc_C(int N, double* x, double *prev_x, double* f, double* g) {
18		✗	GEO::geo_argused(prev_x);
19		✗	*f = 0;
20		✗	for (int i = 0; i < N; i+=2) {
21		✗	double T1 = 1 - x[i];
22		✗	double T2 = 10*(x[i+1]-x[i]*x[i]);
23		✗	*f += T1*T1+T2*T2;
24		✗	g[i+1] = 20*T2;
25		✗	g[i] = -2*(x[i]*g[i+1]+T1);
26			}
27		✗	}
28
29		✗	void newiteration_C(
30			int iter, int call_iter, double *x,
31			double* f, double *g, double* gnorm
32			) {
33		✗	GEO::geo_argused(x);
34		✗	GEO::geo_argused(g);
35		✗	std::cout << iter <<": " << call_iter <<" "
36		✗	<< *f <<" " << *gnorm << std::endl;
37		✗	}
38
39		✗	void evalfunc_h_C(
40			int N, double *x, double *prev_x, double *f, double *g,
41			HESSIAN_MATRIX& hessian
42			) {
43		✗	GEO::geo_argused(prev_x);
44
45			//the following code is not optimal if the pattern of
46			// hessian matrix is fixed.
47		✗	if (m_sparse_matrix) {
48		✗	delete m_sparse_matrix;
49			}
50
51		✗	m_sparse_matrix = new Lite_Sparse_Matrix(
52			(unsigned int)N, (unsigned int)N, SYM_LOWER, CCS, FORTRAN_TYPE, true
53		✗	);
54
55		✗	m_sparse_matrix->begin_fill_entry();
56
57			static bool first = true;
58		✗	double *diag = m_sparse_matrix->get_diag();
59
60		✗	if (first) {
61			// you need to update f and g
62		✗	*f = 0;
63			double tmp;
64		✗	for (unsigned int i = 0; i < (unsigned int)N; i+=2) {
65		✗	tmp = x[i]*x[i];
66		✗	double T1 = 1 - x[i];
67		✗	double T2 = 10*(x[i+1]-tmp);
68		✗	*f += T1*T1+T2*T2;
69		✗	g[i+1] = 20*T2;
70		✗	g[i] = -2*(x[i]*g[i+1]+T1);
71		✗	diag[i] = 2+1200*tmp-400*x[i+1];
72		✗	diag[i+1] = 200;
73		✗	m_sparse_matrix->fill_entry(i, i+1, -400*x[i]);
74			}
75			} else {
76		✗	for (unsigned int i = 0; i < (unsigned int)N; i+=2) {
77		✗	diag[i] = 2+1200*x[i]*x[i]-400*x[i+1];
78		✗	diag[i+1] = 200;
79		✗	m_sparse_matrix->fill_entry(i, i+1, -400*x[i]);
80			}
81			}
82		✗	m_sparse_matrix->end_fill_entry();
83		✗	hessian.set_diag(m_sparse_matrix->get_diag());
84		✗	hessian.set_values(m_sparse_matrix->get_values());
85		✗	hessian.set_rowind((int*)m_sparse_matrix->get_rowind());
86		✗	hessian.set_colptr((int*)m_sparse_matrix->get_colptr());
87		✗	hessian.set_nonzeros((int)m_sparse_matrix->get_nonzero());
88		✗	first = false;
89		✗	}
90
91		✗	void Optimize_by_HLBFGS_C(
92			int N, double *init_x, int num_iter, int M, int T, bool with_hessian
93			) {
94			double parameter[20];
95			int info[20];
96			//initialize
97		✗	INIT_HLBFGS(parameter, info);
98		✗	info[4] = num_iter;
99		✗	info[6] = T;
100		✗	info[7] = with_hessian?1:0;
101		✗	info[10] = 0;
102		✗	info[11] = 1;
103
104		✗	if (with_hessian) {
105		✗	HLBFGS(
106			N, M, init_x,
107			evalfunc_C, evalfunc_h_C,
108			HLBFGS_UPDATE_Hessian, newiteration_C, parameter, info
109			);
110			} else {
111		✗	HLBFGS(
112			N, M, init_x,
113			evalfunc_C, nullptr,
114			HLBFGS_UPDATE_Hessian, newiteration_C, parameter, info
115			);
116			}
117		✗	}
118
119			/********************** Geogram API **************************************/
120
121		✗	void evalfunc(GEO::index_t N, double* x, double& f, double* g) {
122		✗	std::cerr << "eval func" << std::endl;
123		✗	f = 0.0;
124		✗	for (GEO::index_t i = 0; i < N; i+=2) {
125		✗	double T1 = 1 - x[i];
126		✗	double T2 = 10*(x[i+1]-x[i]*x[i]);
127		✗	f += T1*T1+T2*T2;
128		✗	g[i+1] = 20*T2;
129		✗	g[i] = -2*(x[i]*g[i+1]+T1);
130			}
131		✗	}
132
133		✗	void newiteration(
134			GEO::index_t iter,
135			const double *x, double f, const double *g, double gnorm
136			) {
137		✗	GEO::geo_argused(iter);
138		✗	GEO::geo_argused(x);
139		✗	GEO::geo_argused(g);
140		✗	std::cout << " " << f <<" " << gnorm << std::endl;
141		✗	}
142
143		✗	void evalfunc_h(
144			GEO::index_t N, double *x, double& f, double *g,
145			HESSIAN_MATRIX& hessian
146			) {
147		✗	std::cerr << "eval func with Hessian" << std::endl;
148
149			//the following code is not optimal if the pattern of
150			// hessian matrix is fixed.
151		✗	if (m_sparse_matrix) {
152		✗	delete m_sparse_matrix;
153			}
154
155		✗	m_sparse_matrix = new Lite_Sparse_Matrix(
156			N, N, SYM_LOWER, CCS, FORTRAN_TYPE, true
157		✗	);
158
159		✗	m_sparse_matrix->begin_fill_entry();
160
161			static bool first = true;
162		✗	double *diag = m_sparse_matrix->get_diag();
163
164		✗	if (first) {
165			// you need to update f and g
166		✗	f = 0.0;
167			double tmp;
168		✗	for (unsigned int i = 0; i < N; i+=2) {
169		✗	tmp = x[i]*x[i];
170		✗	double T1 = 1 - x[i];
171		✗	double T2 = 10*(x[i+1]-tmp);
172		✗	f += T1*T1+T2*T2;
173		✗	g[i+1] = 20*T2;
174		✗	g[i] = -2*(x[i]*g[i+1]+T1);
175		✗	diag[i] = 2+1200*tmp-400*x[i+1];
176		✗	diag[i+1] = 200;
177		✗	m_sparse_matrix->fill_entry(i, i+1, -400*x[i]);
178			}
179			} else {
180		✗	for (unsigned int i = 0; i < N; i+=2) {
181		✗	diag[i] = 2+1200*x[i]*x[i]-400*x[i+1];
182		✗	diag[i+1] = 200;
183		✗	m_sparse_matrix->fill_entry(i, i+1, -400*x[i]);
184			}
185			}
186		✗	m_sparse_matrix->end_fill_entry();
187		✗	hessian.set_diag(m_sparse_matrix->get_diag());
188		✗	hessian.set_values(m_sparse_matrix->get_values());
189		✗	hessian.set_rowind((int*)m_sparse_matrix->get_rowind());
190		✗	hessian.set_colptr((int*)m_sparse_matrix->get_colptr());
191		✗	hessian.set_nonzeros((int)m_sparse_matrix->get_nonzero());
192		✗	first = false;
193		✗	}
194
195		✗	void Optimize_by_HLBFGS(
196			int N, double *init_x, int num_iter, int M, int T, bool with_hessian
197			) {
198			GEO::Optimizer_var optimizer = with_hessian ?
199		✗	GEO::Optimizer::create("HLBFGS_HESS") :
200		✗	GEO::Optimizer::create("HLBFGS");
201
202
203		✗	optimizer->set_newiteration_callback(newiteration);
204		✗	if(with_hessian) {
205		✗	optimizer->set_evalhessian_callback(evalfunc_h);
206			}
207		✗	optimizer->set_funcgrad_callback(evalfunc);
208
209		✗	optimizer->set_N((unsigned int)N);
210		✗	optimizer->set_M((unsigned int)M);
211		✗	optimizer->set_max_iter((unsigned int)num_iter);
212
213			GEO::HLBFGS_HessOptimizer* hess =
214		✗	dynamic_cast<GEO::HLBFGS_HessOptimizer*>(
215		✗	(GEO::Optimizer*)(optimizer)
216		✗	);
217
218		✗	if(hess != nullptr) {
219		✗	hess->set_T((unsigned int)T);
220			}
221		✗	optimizer->optimize(init_x);
222		✗	}
223			}
224
225			/****************************************************************************/
226
227		✗	int main(int argc, char** argv) {
228
229		✗	GEO::initialize(GEO::GEOGRAM_INSTALL_ALL);
230		✗	GEO::CmdLine::import_arg_group("standard");
231		✗	GEO::CmdLine::declare_arg("Newton",false,"Use Newton solver");
232		✗	GEO::CmdLine::declare_arg("C_api",false,"Use HLBFGS C api");
233		✗	GEO::CmdLine::declare_arg("N", 1000, "Nb variables");
234
235		✗	if(!GEO::CmdLine::parse(argc, argv)) {
236		✗	return 1;
237			}
238
239		✗	std::cout.precision(16);
240		✗	std::cout << std::scientific;
241
242		✗	int N = GEO::CmdLine::get_arg_int("N");
243		✗	std::vector<double> x((unsigned int)N);
244
245		✗	for (unsigned int i = 0; i < (unsigned int)(N/2); i++) {
246		✗	x[2*i] = -1.2;
247		✗	x[2*i+1] = 1.0;
248			}
249
250		✗	int M = 7;
251		✗	int T = 0;
252
253		✗	if(GEO::CmdLine::get_arg_bool("C_api")) {
254		✗	if(GEO::CmdLine::get_arg_bool("Newton")) {
255			//use Hessian
256			// if M = 0, T = 0, it is Newton
257		✗	Optimize_by_HLBFGS_C(N, &x[0], 1000, M, T, true);
258			} else {
259			//without Hessian
260			// it is LBFGS(M) actually, T is not used
261		✗	Optimize_by_HLBFGS_C(N, &x[0], 1000, M, T, false);
262			}
263			} else {
264		✗	if(GEO::CmdLine::get_arg_bool("Newton")) {
265			//use Hessian
266			// if M = 0, T = 0, it is Newton
267		✗	Optimize_by_HLBFGS(N, &x[0], 1000, M, T, true);
268			} else {
269			//without Hessian
270			// it is LBFGS(M) actually, T is not used
271		✗	Optimize_by_HLBFGS(N, &x[0], 1000, M, T, false);
272			}
273			}
274
275		✗	if (m_sparse_matrix) {
276		✗	delete m_sparse_matrix;
277			}
278
279		✗	return 0;
280		✗	}
281
282			#else
283
284			#include <iostream>
285
286			int main() {
287			std::cout << "This geogram was not compiled with HLBFGS support"
288			<< std::endl;
289			return 0;
290			}
291
292			#endif
293