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

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