GCC Code Coverage Report

Line	Branch	Exec	Source
1			/*
2			* Fitting Tools
3			*
4			* Authors:
5			* Marco Cecchetti <mrcekets at gmail.com>
6			*
7			* Copyright 2008 authors
8			*
9			* This library is free software; you can redistribute it and/or
10			* modify it either under the terms of the GNU Lesser General Public
11			* License version 2.1 as published by the Free Software Foundation
12			* (the "LGPL") or, at your option, under the terms of the Mozilla
13			* Public License Version 1.1 (the "MPL"). If you do not alter this
14			* notice, a recipient may use your version of this file under either
15			* the MPL or the LGPL.
16			*
17			* You should have received a copy of the LGPL along with this library
18			* in the file COPYING-LGPL-2.1; if not, write to the Free Software
19			* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20			* You should have received a copy of the MPL along with this library
21			* in the file COPYING-MPL-1.1
22			*
23			* The contents of this file are subject to the Mozilla Public License
24			* Version 1.1 (the "License"); you may not use this file except in
25			* compliance with the License. You may obtain a copy of the License at
26			* http://www.mozilla.org/MPL/
27			*
28			* This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
29			* OF ANY KIND, either express or implied. See the LGPL or the MPL for
30			* the specific language governing rights and limitations.
31			*/
32
33
34			#ifndef _NL_FITTING_TOOL_H_
35			#define _NL_FITTING_TOOL_H_
36
37
38			#include <2geom/numeric/vector.h>
39			#include <2geom/numeric/matrix.h>
40
41			#include <2geom/point.h>
42
43			#include <vector>
44
45
46			/*
47			* The least_square_fitter class represents a tool for solving a fitting
48			* problem with respect to a given model that represents an expression
49			* dependent from a parameter where the coefficients of this expression
50			* are the unknowns of the fitting problem.
51			* The minimizing solution is found by computing the pseudo-inverse
52			* of the problem matrix
53			*/
54
55
56			namespace Geom { namespace NL {
57
58			namespace detail {
59
60			template< typename ModelT>
61			class lsf_base
62			{
63			public:
64			typedef ModelT model_type;
65			typedef typename model_type::parameter_type parameter_type;
66			typedef typename model_type::value_type value_type;
67
68		✗	lsf_base(model_type const &_model, size_t forecasted_samples)
69		✗	: m_model(_model)
70		✗	, m_total_samples(0)
71		✗	, m_matrix(forecasted_samples, m_model.size())
72		✗	, m_psdinv_matrix(nullptr)
73		✗	{}
74
75			// compute pseudo inverse
76		✗	void update()
77			{
78		✗	if (total_samples() == 0) return;
79		✗	if (m_psdinv_matrix != NULL)
80			{
81		✗	delete m_psdinv_matrix;
82			}
83		✗	MatrixView mv(m_matrix, 0, 0, total_samples(), m_matrix.columns());
84		✗	m_psdinv_matrix = new Matrix( pseudo_inverse(mv) );
85		✗	assert(m_psdinv_matrix != NULL);
86		✗	}
87
88		✗	size_t total_samples() const
89			{
90		✗	return m_total_samples;
91			}
92
93		✗	bool is_full() const
94			{
95		✗	return (total_samples() == m_matrix.rows());
96			}
97
98			void clear()
99			{
100			m_total_samples = 0;
101			}
102
103			virtual
104		✗	~lsf_base()
105			{
106		✗	if (m_psdinv_matrix != NULL)
107			{
108		✗	delete m_psdinv_matrix;
109			}
110		✗	}
111
112			protected:
113			const model_type & m_model;
114			size_t m_total_samples;
115			Matrix m_matrix;
116			Matrix* m_psdinv_matrix;
117
118			}; // end class lsf_base
119
120
121
122
123			template< typename ModelT, typename ValueType = typename ModelT::value_type>
124			class lsf_solution
125			{
126			};
127
128			// a fitting process on samples with value of type double
129			// produces a solution of type Vector
130			template< typename ModelT>
131			class lsf_solution<ModelT, double>
132			: public lsf_base<ModelT>
133			{
134			public:
135			typedef ModelT model_type;
136			typedef typename model_type::parameter_type parameter_type;
137			typedef typename model_type::value_type value_type;
138			typedef Vector solution_type;
139			typedef lsf_base<model_type> base_type;
140
141			using base_type::m_model;
142			using base_type::m_psdinv_matrix;
143			using base_type::total_samples;
144
145			public:
146		✗	lsf_solution(model_type const &_model,
147			size_t forecasted_samples)
148			: base_type(_model, forecasted_samples)
149		✗	, m_solution(_model.size())
150		✗	{}
151
152			template< typename VectorT >
153		✗	solution_type& result(VectorT const& sample_values)
154			{
155		✗	assert(sample_values.size() == total_samples());
156		✗	ConstVectorView sv(sample_values);
157		✗	m_solution = (*m_psdinv_matrix) * sv;
158		✗	return m_solution;
159		✗	}
160
161			// a comparison between old sample values and the new ones is performed
162			// in order to minimize computation
163			// prerequisite:
164			// old_sample_values.size() == new_sample_values.size()
165			// no update() call can be performed between two result invocations
166			template< typename VectorT >
167			solution_type& result( VectorT const& old_sample_values,
168			VectorT const& new_sample_values )
169			{
170			assert(old_sample_values.size() == total_samples());
171			assert(new_sample_values.size() == total_samples());
172			Vector diff(total_samples());
173			for (size_t i = 0; i < diff.size(); ++i)
174			{
175			diff[i] = new_sample_values[i] - old_sample_values[i];
176			}
177			Vector column(m_model.size());
178			Vector delta(m_model.size(), 0.0);
179			for (size_t i = 0; i < diff.size(); ++i)
180			{
181			if (diff[i] != 0)
182			{
183			column = m_psdinv_matrix->column_view(i);
184			column.scale(diff[i]);
185			delta += column;
186			}
187			}
188			m_solution += delta;
189			return m_solution;
190			}
191
192		✗	solution_type& result()
193			{
194		✗	return m_solution;
195			}
196
197			private:
198			solution_type m_solution;
199
200			}; // end class lsf_solution<ModelT, double>
201
202
203			// a fitting process on samples with value of type Point
204			// produces a solution of type Matrix (with 2 columns)
205			template< typename ModelT>
206			class lsf_solution<ModelT, Point>
207			: public lsf_base<ModelT>
208			{
209			public:
210			typedef ModelT model_type;
211			typedef typename model_type::parameter_type parameter_type;
212			typedef typename model_type::value_type value_type;
213			typedef Matrix solution_type;
214			typedef lsf_base<model_type> base_type;
215
216			using base_type::m_model;
217			using base_type::m_psdinv_matrix;
218			using base_type::total_samples;
219
220			public:
221			lsf_solution(model_type const &_model,
222			size_t forecasted_samples)
223			: base_type(_model, forecasted_samples)
224			, m_solution(_model.size(), 2)
225			{}
226
227			solution_type& result(std::vector<Point> const& sample_values)
228			{
229			assert(sample_values.size() == total_samples());
230			Matrix svm(total_samples(), 2);
231			for (size_t i = 0; i < total_samples(); ++i)
232			{
233			svm(i, X) = sample_values[i][X];
234			svm(i, Y) = sample_values[i][Y];
235			}
236			m_solution = (*m_psdinv_matrix) * svm;
237			return m_solution;
238			}
239
240			// a comparison between old sample values and the new ones is performed
241			// in order to minimize computation
242			// prerequisite:
243			// old_sample_values.size() == new_sample_values.size()
244			// no update() call can to be performed between two result invocations
245			solution_type& result( std::vector<Point> const& old_sample_values,
246			std::vector<Point> const& new_sample_values )
247			{
248			assert(old_sample_values.size() == total_samples());
249			assert(new_sample_values.size() == total_samples());
250			Matrix diff(total_samples(), 2);
251			for (size_t i = 0; i < total_samples(); ++i)
252			{
253			diff(i, X) = new_sample_values[i][X] - old_sample_values[i][X];
254			diff(i, Y) = new_sample_values[i][Y] - old_sample_values[i][Y];
255			}
256			Vector column(m_model.size());
257			Matrix delta(m_model.size(), 2, 0.0);
258			VectorView deltax = delta.column_view(X);
259			VectorView deltay = delta.column_view(Y);
260			for (size_t i = 0; i < total_samples(); ++i)
261			{
262			if (diff(i, X) != 0)
263			{
264			column = m_psdinv_matrix->column_view(i);
265			column.scale(diff(i, X));
266			deltax += column;
267			}
268			if (diff(i, Y) != 0)
269			{
270			column = m_psdinv_matrix->column_view(i);
271			column.scale(diff(i, Y));
272			deltay += column;
273			}
274			}
275			m_solution += delta;
276			return m_solution;
277			}
278
279			solution_type& result()
280			{
281			return m_solution;
282			}
283
284			private:
285			solution_type m_solution;
286
287			}; // end class lsf_solution<ModelT, Point>
288
289
290
291
292			template< typename ModelT,
293			bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
294			class lsf_with_fixed_terms
295			{
296			};
297
298
299			// fitting tool for completely unknown models
300			template< typename ModelT>
301			class lsf_with_fixed_terms<ModelT, false>
302			: public lsf_solution<ModelT>
303			{
304			public:
305			typedef ModelT model_type;
306			typedef typename model_type::parameter_type parameter_type;
307			typedef typename model_type::value_type value_type;
308			typedef lsf_solution<model_type> base_type;
309			typedef typename base_type::solution_type solution_type;
310
311			using base_type::total_samples;
312			using base_type::is_full;
313			using base_type::m_matrix;
314			using base_type::m_total_samples;
315			using base_type::m_model;
316
317			public:
318			lsf_with_fixed_terms(model_type const &_model,
319			size_t forecasted_samples)
320			: base_type(_model, forecasted_samples)
321			{}
322
323			void append(parameter_type const& sample_parameter)
324			{
325			assert(!is_full());
326			VectorView row = m_matrix.row_view(total_samples());
327			m_model.feed(row, sample_parameter);
328			++m_total_samples;
329			}
330
331			void append_copy(size_t sample_index)
332			{
333			assert(!is_full());
334			assert(sample_index < total_samples());
335			VectorView dest_row = m_matrix.row_view(total_samples());
336			VectorView source_row = m_matrix.row_view(sample_index);
337			dest_row = source_row;
338			++m_total_samples;
339			}
340
341			}; // end class lsf_with_fixed_terms<ModelT, false>
342
343
344			// fitting tool for partially known models
345			template< typename ModelT>
346			class lsf_with_fixed_terms<ModelT, true>
347			: public lsf_solution<ModelT>
348			{
349			public:
350			typedef ModelT model_type;
351			typedef typename model_type::parameter_type parameter_type;
352			typedef typename model_type::value_type value_type;
353			typedef lsf_solution<model_type> base_type;
354			typedef typename base_type::solution_type solution_type;
355
356			using base_type::total_samples;
357			using base_type::is_full;
358			using base_type::m_matrix;
359			using base_type::m_total_samples;
360			using base_type::m_model;
361
362			public:
363		✗	lsf_with_fixed_terms(model_type const &_model,
364			size_t forecasted_samples)
365			: base_type(_model, forecasted_samples)
366		✗	, m_vector(forecasted_samples)
367		✗	, m_vector_view(nullptr)
368		✗	{}
369
370		✗	void append(parameter_type const& sample_parameter)
371			{
372		✗	assert(!is_full());
373		✗	VectorView row = m_matrix.row_view(total_samples());
374		✗	m_model.feed(row, m_vector[total_samples()], sample_parameter);
375		✗	++m_total_samples;
376		✗	}
377
378			void append_copy(size_t sample_index)
379			{
380			assert(!is_full());
381			assert(sample_index < total_samples());
382			VectorView dest_row = m_matrix.row_view(total_samples());
383			VectorView source_row = m_matrix.row_view(sample_index);
384			dest_row = source_row;
385			m_vector[total_samples()] = m_vector[sample_index];
386			++m_total_samples;
387			}
388
389		✗	void update()
390			{
391		✗	base_type::update();
392		✗	if (total_samples() == 0) return;
393		✗	if (m_vector_view != NULL)
394			{
395		✗	delete m_vector_view;
396			}
397		✗	m_vector_view = new VectorView(m_vector, base_type::total_samples());
398		✗	assert(m_vector_view != NULL);
399			}
400
401
402		✗	~lsf_with_fixed_terms() override
403			{
404		✗	if (m_vector_view) {
405		✗	delete m_vector_view;
406			}
407		✗	}
408
409			protected:
410			Vector m_vector;
411			VectorView* m_vector_view;
412			};
413
414
415			} // end namespace detail
416
417
418
419
420			template< typename ModelT,
421			typename ValueType = typename ModelT::value_type,
422			bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
423			class least_squeares_fitter
424			{
425			};
426
427
428			template< typename ModelT, typename ValueType >
429			class least_squeares_fitter<ModelT, ValueType, false>
430			: public detail::lsf_with_fixed_terms<ModelT>
431			{
432			public:
433			typedef ModelT model_type;
434			typedef detail::lsf_with_fixed_terms<model_type> base_type;
435			typedef typename base_type::parameter_type parameter_type;
436			typedef typename base_type::value_type value_type;
437			typedef typename base_type::solution_type solution_type;
438
439			public:
440			least_squeares_fitter(model_type const &_model,
441			size_t forecasted_samples)
442			: base_type(_model, forecasted_samples)
443			{}
444			};
445
446			template< typename ModelT>
447			class least_squeares_fitter<ModelT, double, true>
448			: public detail::lsf_with_fixed_terms<ModelT>
449			{
450			public:
451			typedef ModelT model_type;
452			typedef detail::lsf_with_fixed_terms<model_type> base_type;
453			typedef typename base_type::parameter_type parameter_type;
454			typedef typename base_type::value_type value_type;
455			typedef typename base_type::solution_type solution_type;
456
457			using base_type::m_vector_view;
458			//using base_type::result; // VSC legacy support
459			solution_type& result( std::vector<Point> const& old_sample_values,
460			std::vector<Point> const& new_sample_values )
461			{
462			return base_type::result(old_sample_values, new_sample_values);
463			}
464
465		✗	solution_type& result()
466			{
467		✗	return base_type::result();
468			}
469
470			public:
471		✗	least_squeares_fitter(model_type const &_model,
472			size_t forecasted_samples)
473		✗	: base_type(_model, forecasted_samples)
474		✗	{}
475
476			template< typename VectorT >
477		✗	solution_type& result(VectorT const& sample_values)
478			{
479		✗	assert(sample_values.size() == m_vector_view->size());
480		✗	Vector sv(sample_values.size());
481		✗	for (size_t i = 0; i < sv.size(); ++i)
482		✗	sv[i] = sample_values[i] - (*m_vector_view)[i];
483		✗	return base_type::result(sv);
484		✗	}
485
486			}; // end class least_squeares_fitter<ModelT, double, true>
487
488
489			template< typename ModelT>
490			class least_squeares_fitter<ModelT, Point, true>
491			: public detail::lsf_with_fixed_terms<ModelT>
492			{
493			public:
494			typedef ModelT model_type;
495			typedef detail::lsf_with_fixed_terms<model_type> base_type;
496			typedef typename base_type::parameter_type parameter_type;
497			typedef typename base_type::value_type value_type;
498			typedef typename base_type::solution_type solution_type;
499
500			using base_type::m_vector_view;
501			//using base_type::result; // VCS legacy support
502			solution_type& result( std::vector<Point> const& old_sample_values,
503			std::vector<Point> const& new_sample_values )
504			{
505			return base_type::result(old_sample_values, new_sample_values);
506			}
507
508			solution_type& result()
509			{
510			return base_type::result();
511			}
512
513
514			public:
515			least_squeares_fitter(model_type const &_model,
516			size_t forecasted_samples)
517			: base_type(_model, forecasted_samples)
518			{}
519
520			solution_type& result(std::vector<Point> const& sample_values)
521			{
522			assert(sample_values.size() == m_vector_view->size());
523			NL::Matrix sv(sample_values.size(), 2);
524			for (size_t i = 0; i < sample_values.size(); ++i)
525			{
526			sv(i, X) = sample_values[i][X] - (*m_vector_view)[i];
527			sv(i, Y) = sample_values[i][Y] - (*m_vector_view)[i];
528			}
529			return base_type::result(sv);
530			}
531
532			}; // end class least_squeares_fitter<ModelT, Point, true>
533
534
535			} // end namespace NL
536			} // end namespace Geom
537
538
539
540			#endif // _NL_FITTING_TOOL_H_
541
542
543			/*
544			Local Variables:
545			mode:c++
546			c-file-style:"stroustrup"
547			c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
548			indent-tabs-mode:nil
549			fill-column:99
550			End:
551			*/
552			// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
553