src/2geom/numeric/fitting-tool.h

   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(NULL)
  73     {
  74     }
  75
  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     }
  88
  89     size_t total_samples() const
  90     {
  91         return m_total_samples;
  92     }
  93
  94     bool is_full() const
  95     {
  96         return (total_samples() == m_matrix.rows());
  97     }
  98
  99     void clear()
 100     {
 101         m_total_samples = 0;
 102     }
 103
 104     virtual
 105     ~lsf_base()
 106     {
 107         if (m_psdinv_matrix != NULL)
 108         {
 109             delete m_psdinv_matrix;
 110         }
 111     }
 112
 113   protected:
 114     const model_type &      m_model;
 115     size_t                  m_total_samples;
 116     Matrix                  m_matrix;
 117     Matrix*                 m_psdinv_matrix;
 118
 119 };  // end class lsf_base
 120
 121
 122
 123
 124 template< typename ModelT, typename ValueType = typename ModelT::value_type>
 125 class lsf_solution
 126 {
 127 };
 128
 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;
 141
 142     using base_type::m_model;
 143     using base_type::m_psdinv_matrix;
 144     using base_type::total_samples;
 145
 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     }
 153
 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     }
 162
 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     }
 193
 194     solution_type& result()
 195     {
 196         return m_solution;
 197     }
 198
 199 private:
 200     solution_type m_solution;
 201
 202 }; // end class lsf_solution<ModelT, double>
 203
 204
 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;
 217
 218     using base_type::m_model;
 219     using base_type::m_psdinv_matrix;
 220     using base_type::total_samples;
 221
 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     }
 229
 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     }
 242
 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     }
 281
 282     solution_type& result()
 283     {
 284         return m_solution;
 285     }
 286
 287 private:
 288     solution_type m_solution;
 289
 290 }; // end class lsf_solution<ModelT, Point>
 291
 292
 293
 294
 295 template< typename ModelT,
 296           bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
 297 class lsf_with_fixed_terms
 298 {
 299 };
 300
 301
 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;
 313
 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;
 319
 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     }
 326
 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     }
 334
 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     }
 344
 345 }; // end class lsf_with_fixed_terms<ModelT, false>
 346
 347
 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;
 359
 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;
 365
 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     }
 381
 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     }
 392
 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     }
 404
 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     }
 413
 414   protected:
 415     Vector          m_vector;
 416     VectorView*     m_vector_view;
 417
 418 }; // end class lsf_with_fixed_terms<ModelT, true>
 419
 420
 421 } // end namespace detail
 422
 423
 424
 425
 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 };
 432
 433
 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;
 444
 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>
 452
 453
 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;
 464
 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     }
 472
 473     solution_type& result()
 474     {
 475         return base_type::result();
 476     }
 477
 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     }
 484
 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     }
 494
 495 }; // end class least_squeares_fitter<ModelT, double, true>
 496
 497
 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;
 508
 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     }
 516
 517     solution_type& result()
 518     {
 519         return base_type::result();
 520     }
 521
 522
 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     }
 529
 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     }
 541
 542 };  // end class least_squeares_fitter<ModelT, Point, true>
 543
 544
 545 }  // end namespace NL
 546 }  // end namespace Geom
 547
 548
 549
 550 #endif  // _NL_FITTING_TOOL_H_
 551
 552
 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 :