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 */
34 #ifndef _NL_FITTING_TOOL_H_
35 #define _NL_FITTING_TOOL_H_
38 #include <2geom/numeric/vector.h>
39 #include <2geom/numeric/matrix.h>
41 #include <2geom/point.h>
43 #include <vector>
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 */
56 namespace Geom { namespace NL {
58 namespace detail {
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;
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(NULL)
73 {
74 }
76 // compute pseudo inverse
77 void update()
78 {
79 if (total_samples() == 0) return;
80 if (m_psdinv_matrix != NULL)
81 {
82 delete m_psdinv_matrix;
83 }
84 MatrixView mv(m_matrix, 0, 0, total_samples(), m_matrix.columns());
85 m_psdinv_matrix = new Matrix( pseudo_inverse(mv) );
86 assert(m_psdinv_matrix != NULL);
87 }
89 size_t total_samples() const
90 {
91 return m_total_samples;
92 }
94 bool is_full() const
95 {
96 return (total_samples() == m_matrix.rows());
97 }
99 void clear()
100 {
101 m_total_samples = 0;
102 }
104 virtual
105 ~lsf_base()
106 {
107 if (m_psdinv_matrix != NULL)
108 {
109 delete m_psdinv_matrix;
110 }
111 }
113 protected:
114 const model_type & m_model;
115 size_t m_total_samples;
116 Matrix m_matrix;
117 Matrix* m_psdinv_matrix;
119 }; // end class lsf_base
124 template< typename ModelT, typename ValueType = typename ModelT::value_type>
125 class lsf_solution
126 {
127 };
129 // a fitting process on samples with value of type double
130 // produces a solution of type Vector
131 template< typename ModelT>
132 class lsf_solution<ModelT, double>
133 : public lsf_base<ModelT>
134 {
135 public:
136 typedef ModelT model_type;
137 typedef typename model_type::parameter_type parameter_type;
138 typedef typename model_type::value_type value_type;
139 typedef Vector solution_type;
140 typedef lsf_base<model_type> base_type;
142 using base_type::m_model;
143 using base_type::m_psdinv_matrix;
144 using base_type::total_samples;
146 public:
147 lsf_solution<ModelT, double>( model_type const& _model,
148 size_t forecasted_samples )
149 : base_type(_model, forecasted_samples),
150 m_solution(_model.size())
151 {
152 }
154 template< typename VectorT >
155 solution_type& result(VectorT const& sample_values)
156 {
157 assert(sample_values.size() == total_samples());
158 ConstVectorView sv(sample_values);
159 m_solution = (*m_psdinv_matrix) * sv;
160 return m_solution;
161 }
163 // a comparison between old sample values and the new ones is performed
164 // in order to minimize computation
165 // prerequisite:
166 // old_sample_values.size() == new_sample_values.size()
167 // no update() call can be performed between two result invocations
168 template< typename VectorT >
169 solution_type& result( VectorT const& old_sample_values,
170 VectorT const& new_sample_values )
171 {
172 assert(old_sample_values.size() == total_samples());
173 assert(new_sample_values.size() == total_samples());
174 Vector diff(total_samples());
175 for (size_t i = 0; i < diff.size(); ++i)
176 {
177 diff[i] = new_sample_values[i] - old_sample_values[i];
178 }
179 Vector column(m_model.size());
180 Vector delta(m_model.size(), 0.0);
181 for (size_t i = 0; i < diff.size(); ++i)
182 {
183 if (diff[i] != 0)
184 {
185 column = m_psdinv_matrix->column_view(i);
186 column.scale(diff[i]);
187 delta += column;
188 }
189 }
190 m_solution += delta;
191 return m_solution;
192 }
194 solution_type& result()
195 {
196 return m_solution;
197 }
199 private:
200 solution_type m_solution;
202 }; // end class lsf_solution<ModelT, double>
205 // a fitting process on samples with value of type Point
206 // produces a solution of type Matrix (with 2 columns)
207 template< typename ModelT>
208 class lsf_solution<ModelT, Point>
209 : public lsf_base<ModelT>
210 {
211 public:
212 typedef ModelT model_type;
213 typedef typename model_type::parameter_type parameter_type;
214 typedef typename model_type::value_type value_type;
215 typedef Matrix solution_type;
216 typedef lsf_base<model_type> base_type;
218 using base_type::m_model;
219 using base_type::m_psdinv_matrix;
220 using base_type::total_samples;
222 public:
223 lsf_solution<ModelT, Point>( model_type const& _model,
224 size_t forecasted_samples )
225 : base_type(_model, forecasted_samples),
226 m_solution(_model.size(), 2)
227 {
228 }
230 solution_type& result(std::vector<Point> const& sample_values)
231 {
232 assert(sample_values.size() == total_samples());
233 Matrix svm(total_samples(), 2);
234 for (size_t i = 0; i < total_samples(); ++i)
235 {
236 svm(i, X) = sample_values[i][X];
237 svm(i, Y) = sample_values[i][Y];
238 }
239 m_solution = (*m_psdinv_matrix) * svm;
240 return m_solution;
241 }
243 // a comparison between old sample values and the new ones is performed
244 // in order to minimize computation
245 // prerequisite:
246 // old_sample_values.size() == new_sample_values.size()
247 // no update() call can to be performed between two result invocations
248 solution_type& result( std::vector<Point> const& old_sample_values,
249 std::vector<Point> const& new_sample_values )
250 {
251 assert(old_sample_values.size() == total_samples());
252 assert(new_sample_values.size() == total_samples());
253 Matrix diff(total_samples(), 2);
254 for (size_t i = 0; i < total_samples(); ++i)
255 {
256 diff(i, X) = new_sample_values[i][X] - old_sample_values[i][X];
257 diff(i, Y) = new_sample_values[i][Y] - old_sample_values[i][Y];
258 }
259 Vector column(m_model.size());
260 Matrix delta(m_model.size(), 2, 0.0);
261 VectorView deltax = delta.column_view(X);
262 VectorView deltay = delta.column_view(Y);
263 for (size_t i = 0; i < total_samples(); ++i)
264 {
265 if (diff(i, X) != 0)
266 {
267 column = m_psdinv_matrix->column_view(i);
268 column.scale(diff(i, X));
269 deltax += column;
270 }
271 if (diff(i, Y) != 0)
272 {
273 column = m_psdinv_matrix->column_view(i);
274 column.scale(diff(i, Y));
275 deltay += column;
276 }
277 }
278 m_solution += delta;
279 return m_solution;
280 }
282 solution_type& result()
283 {
284 return m_solution;
285 }
287 private:
288 solution_type m_solution;
290 }; // end class lsf_solution<ModelT, Point>
295 template< typename ModelT,
296 bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
297 class lsf_with_fixed_terms
298 {
299 };
302 // fitting tool for completely unknown models
303 template< typename ModelT>
304 class lsf_with_fixed_terms<ModelT, false>
305 : public lsf_solution<ModelT>
306 {
307 public:
308 typedef ModelT model_type;
309 typedef typename model_type::parameter_type parameter_type;
310 typedef typename model_type::value_type value_type;
311 typedef lsf_solution<model_type> base_type;
312 typedef typename base_type::solution_type solution_type;
314 using base_type::total_samples;
315 using base_type::is_full;
316 using base_type::m_matrix;
317 using base_type::m_total_samples;
318 using base_type::m_model;
320 public:
321 lsf_with_fixed_terms<ModelT, false>( model_type const& _model,
322 size_t forecasted_samples )
323 : base_type(_model, forecasted_samples)
324 {
325 }
327 void append(parameter_type const& sample_parameter)
328 {
329 assert(!is_full());
330 VectorView row = m_matrix.row_view(total_samples());
331 m_model.feed(row, sample_parameter);
332 ++m_total_samples;
333 }
335 void append_copy(size_t sample_index)
336 {
337 assert(!is_full());
338 assert(sample_index < total_samples());
339 VectorView dest_row = m_matrix.row_view(total_samples());
340 VectorView source_row = m_matrix.row_view(sample_index);
341 dest_row = source_row;
342 ++m_total_samples;
343 }
345 }; // end class lsf_with_fixed_terms<ModelT, false>
348 // fitting tool for partially known models
349 template< typename ModelT>
350 class lsf_with_fixed_terms<ModelT, true>
351 : public lsf_solution<ModelT>
352 {
353 public:
354 typedef ModelT model_type;
355 typedef typename model_type::parameter_type parameter_type;
356 typedef typename model_type::value_type value_type;
357 typedef lsf_solution<model_type> base_type;
358 typedef typename base_type::solution_type solution_type;
360 using base_type::total_samples;
361 using base_type::is_full;
362 using base_type::m_matrix;
363 using base_type::m_total_samples;
364 using base_type::m_model;
366 public:
367 lsf_with_fixed_terms<ModelT, true>( model_type const& _model,
368 size_t forecasted_samples )
369 : base_type(_model, forecasted_samples),
370 m_vector(forecasted_samples),
371 m_vector_view(NULL)
372 {
373 }
374 void append(parameter_type const& sample_parameter)
375 {
376 assert(!is_full());
377 VectorView row = m_matrix.row_view(total_samples());
378 m_model.feed(row, m_vector[total_samples()], sample_parameter);
379 ++m_total_samples;
380 }
382 void append_copy(size_t sample_index)
383 {
384 assert(!is_full());
385 assert(sample_index < total_samples());
386 VectorView dest_row = m_matrix.row_view(total_samples());
387 VectorView source_row = m_matrix.row_view(sample_index);
388 dest_row = source_row;
389 m_vector[total_samples()] = m_vector[sample_index];
390 ++m_total_samples;
391 }
393 void update()
394 {
395 base_type::update();
396 if (total_samples() == 0) return;
397 if (m_vector_view != NULL)
398 {
399 delete m_vector_view;
400 }
401 m_vector_view = new VectorView(m_vector, base_type::total_samples());
402 assert(m_vector_view != NULL);
403 }
405 virtual
406 ~lsf_with_fixed_terms<model_type, true>()
407 {
408 if (m_vector_view != NULL)
409 {
410 delete m_vector_view;
411 }
412 }
414 protected:
415 Vector m_vector;
416 VectorView* m_vector_view;
418 }; // end class lsf_with_fixed_terms<ModelT, true>
421 } // end namespace detail
426 template< typename ModelT,
427 typename ValueType = typename ModelT::value_type,
428 bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
429 class least_squeares_fitter
430 {
431 };
434 template< typename ModelT, typename ValueType >
435 class least_squeares_fitter<ModelT, ValueType, false>
436 : public detail::lsf_with_fixed_terms<ModelT>
437 {
438 public:
439 typedef ModelT model_type;
440 typedef detail::lsf_with_fixed_terms<model_type> base_type;
441 typedef typename base_type::parameter_type parameter_type;
442 typedef typename base_type::value_type value_type;
443 typedef typename base_type::solution_type solution_type;
445 public:
446 least_squeares_fitter<ModelT, ValueType, false>( model_type const& _model,
447 size_t forecasted_samples )
448 : base_type(_model, forecasted_samples)
449 {
450 }
451 }; // end class least_squeares_fitter<ModelT, ValueType, true>
454 template< typename ModelT>
455 class least_squeares_fitter<ModelT, double, true>
456 : public detail::lsf_with_fixed_terms<ModelT>
457 {
458 public:
459 typedef ModelT model_type;
460 typedef detail::lsf_with_fixed_terms<model_type> base_type;
461 typedef typename base_type::parameter_type parameter_type;
462 typedef typename base_type::value_type value_type;
463 typedef typename base_type::solution_type solution_type;
465 using base_type::m_vector_view;
466 //using base_type::result; // VSC legacy support
467 solution_type& result( std::vector<Point> const& old_sample_values,
468 std::vector<Point> const& new_sample_values )
469 {
470 return base_type::result(old_sample_values, new_sample_values);
471 }
473 solution_type& result()
474 {
475 return base_type::result();
476 }
478 public:
479 least_squeares_fitter<ModelT, double, true>( model_type const& _model,
480 size_t forecasted_samples )
481 : base_type(_model, forecasted_samples)
482 {
483 }
485 template< typename VectorT >
486 solution_type& result(VectorT const& sample_values)
487 {
488 assert(sample_values.size() == m_vector_view->size());
489 Vector sv(sample_values.size());
490 for (size_t i = 0; i < sv.size(); ++i)
491 sv[i] = sample_values[i] - (*m_vector_view)[i];
492 return base_type::result(sv);
493 }
495 }; // end class least_squeares_fitter<ModelT, double, true>
498 template< typename ModelT>
499 class least_squeares_fitter<ModelT, Point, true>
500 : public detail::lsf_with_fixed_terms<ModelT>
501 {
502 public:
503 typedef ModelT model_type;
504 typedef detail::lsf_with_fixed_terms<model_type> base_type;
505 typedef typename base_type::parameter_type parameter_type;
506 typedef typename base_type::value_type value_type;
507 typedef typename base_type::solution_type solution_type;
509 using base_type::m_vector_view;
510 //using base_type::result; // VCS legacy support
511 solution_type& result( std::vector<Point> const& old_sample_values,
512 std::vector<Point> const& new_sample_values )
513 {
514 return base_type::result(old_sample_values, new_sample_values);
515 }
517 solution_type& result()
518 {
519 return base_type::result();
520 }
523 public:
524 least_squeares_fitter<ModelT, Point, true>( model_type const& _model,
525 size_t forecasted_samples )
526 : base_type(_model, forecasted_samples)
527 {
528 }
530 solution_type& result(std::vector<Point> const& sample_values)
531 {
532 assert(sample_values.size() == m_vector_view->size());
533 NL::Matrix sv(sample_values.size(), 2);
534 for (size_t i = 0; i < sample_values.size(); ++i)
535 {
536 sv(i, X) = sample_values[i][X] - (*m_vector_view)[i];
537 sv(i, Y) = sample_values[i][Y] - (*m_vector_view)[i];
538 }
539 return base_type::result(sv);
540 }
542 }; // end class least_squeares_fitter<ModelT, Point, true>
545 } // end namespace NL
546 } // end namespace Geom
550 #endif // _NL_FITTING_TOOL_H_
553 /*
554 Local Variables:
555 mode:c++
556 c-file-style:"stroustrup"
557 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
558 indent-tabs-mode:nil
559 fill-column:99
560 End:
561 */
562 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :