Adding upstream version 1.3.upstream/1.3upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

1 files changed, 562 insertions, 0 deletions

diff --git a/include/2geom/numeric/fitting-tool.h b/include/2geom/numeric/fitting-tool.h
new file mode 100644
index 0000000..78d66ca
--- /dev/null
+++ b/include/2geom/numeric/fitting-tool.h
@@ -0,0 +1,562 @@
+/*
+ * Fitting Tools
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _NL_FITTING_TOOL_H_
+#define _NL_FITTING_TOOL_H_
+
+
+#include <2geom/numeric/vector.h>
+#include <2geom/numeric/matrix.h>
+
+#include <2geom/point.h>
+
+#include <vector>
+
+
+/*
+ *  The least_square_fitter class represents a tool for solving a fitting
+ *  problem with respect to a given model that represents an expression
+ *  dependent from a parameter where the coefficients of this expression
+ *  are the unknowns of the fitting problem.
+ *  The minimizing solution is found by computing the pseudo-inverse
+ *  of the problem matrix
+ */
+
+
+namespace Geom { namespace NL {
+
+namespace detail {
+
+template< typename ModelT>
+class lsf_base
+{
+  public:
+    typedef ModelT                                   model_type;
+    typedef typename model_type::parameter_type      parameter_type;
+    typedef typename model_type::value_type          value_type;
+
+    lsf_base( model_type const& _model, size_t forecasted_samples )
+        : m_model(_model),
+          m_total_samples(0),
+          m_matrix(forecasted_samples, m_model.size()),
+          m_psdinv_matrix(NULL)
+    {
+    }
+
+    // compute pseudo inverse
+    void update()
+    {
+        if (total_samples() == 0) return;
+        if (m_psdinv_matrix != NULL)
+        {
+            delete m_psdinv_matrix;
+        }
+        MatrixView mv(m_matrix, 0, 0, total_samples(), m_matrix.columns());
+        m_psdinv_matrix = new Matrix( pseudo_inverse(mv) );
+        assert(m_psdinv_matrix != NULL);
+    }
+
+    size_t total_samples() const
+    {
+        return m_total_samples;
+    }
+
+    bool is_full() const
+    {
+        return (total_samples() == m_matrix.rows());
+    }
+
+    void clear()
+    {
+        m_total_samples = 0;
+    }
+
+    virtual
+    ~lsf_base()
+    {
+        if (m_psdinv_matrix != NULL)
+        {
+            delete m_psdinv_matrix;
+        }
+    }
+
+  protected:
+    const model_type &      m_model;
+    size_t                  m_total_samples;
+    Matrix                  m_matrix;
+    Matrix*                 m_psdinv_matrix;
+
+};  // end class lsf_base
+
+
+
+
+template< typename ModelT, typename ValueType = typename ModelT::value_type>
+class lsf_solution
+{
+};
+
+// a fitting process on samples with value of type double
+// produces a solution of type Vector
+template< typename ModelT>
+class lsf_solution<ModelT, double>
+    : public lsf_base<ModelT>
+{
+public:
+    typedef ModelT                                  model_type;
+    typedef typename model_type::parameter_type     parameter_type;
+    typedef typename model_type::value_type         value_type;
+    typedef Vector                                  solution_type;
+    typedef lsf_base<model_type>                    base_type;
+
+    using base_type::m_model;
+    using base_type::m_psdinv_matrix;
+    using base_type::total_samples;
+
+public:
+    lsf_solution<ModelT, double>( model_type const& _model,
+                                  size_t forecasted_samples )
+        : base_type(_model, forecasted_samples),
+          m_solution(_model.size())
+    {
+    }
+
+    template< typename VectorT >
+    solution_type& result(VectorT const& sample_values)
+    {
+        assert(sample_values.size() == total_samples());
+        ConstVectorView sv(sample_values);
+        m_solution = (*m_psdinv_matrix) * sv;
+        return m_solution;
+    }
+
+    // a comparison between old sample values and the new ones is performed
+    // in order to minimize computation
+    // prerequisite:
+    //     old_sample_values.size() == new_sample_values.size()
+    //     no update() call can be performed between two result invocations
+    template< typename VectorT >
+    solution_type& result( VectorT const& old_sample_values,
+                           VectorT const& new_sample_values )
+    {
+        assert(old_sample_values.size() == total_samples());
+        assert(new_sample_values.size() == total_samples());
+        Vector diff(total_samples());
+        for (size_t i = 0; i < diff.size(); ++i)
+        {
+            diff[i] = new_sample_values[i] - old_sample_values[i];
+        }
+        Vector column(m_model.size());
+        Vector delta(m_model.size(), 0.0);
+        for (size_t i = 0; i < diff.size(); ++i)
+        {
+            if (diff[i] != 0)
+            {
+                column = m_psdinv_matrix->column_view(i);
+                column.scale(diff[i]);
+                delta += column;
+            }
+        }
+        m_solution += delta;
+        return m_solution;
+    }
+
+    solution_type& result()
+    {
+        return m_solution;
+    }
+
+private:
+    solution_type m_solution;
+
+}; // end class lsf_solution<ModelT, double>
+
+
+// a fitting process on samples with value of type Point
+// produces a solution of type Matrix (with 2 columns)
+template< typename ModelT>
+class lsf_solution<ModelT, Point>
+    : public lsf_base<ModelT>
+{
+public:
+    typedef ModelT                                  model_type;
+    typedef typename model_type::parameter_type     parameter_type;
+    typedef typename model_type::value_type         value_type;
+    typedef Matrix                                  solution_type;
+    typedef lsf_base<model_type>                    base_type;
+
+    using base_type::m_model;
+    using base_type::m_psdinv_matrix;
+    using base_type::total_samples;
+
+public:
+    lsf_solution<ModelT, Point>( model_type const& _model,
+                                 size_t forecasted_samples )
+        : base_type(_model, forecasted_samples),
+          m_solution(_model.size(), 2)
+    {
+    }
+
+    solution_type& result(std::vector<Point> const& sample_values)
+    {
+        assert(sample_values.size() == total_samples());
+        Matrix svm(total_samples(), 2);
+        for (size_t i = 0; i < total_samples(); ++i)
+        {
+            svm(i, X) = sample_values[i][X];
+            svm(i, Y) = sample_values[i][Y];
+        }
+        m_solution = (*m_psdinv_matrix) * svm;
+        return m_solution;
+    }
+
+    // a comparison between old sample values and the new ones is performed
+    // in order to minimize computation
+    // prerequisite:
+    //     old_sample_values.size() == new_sample_values.size()
+    //     no update() call can to be performed between two result invocations
+    solution_type& result( std::vector<Point> const& old_sample_values,
+                           std::vector<Point> const& new_sample_values )
+    {
+        assert(old_sample_values.size() == total_samples());
+        assert(new_sample_values.size() == total_samples());
+        Matrix diff(total_samples(), 2);
+        for (size_t i = 0; i < total_samples(); ++i)
+        {
+            diff(i, X) =  new_sample_values[i][X] - old_sample_values[i][X];
+            diff(i, Y) =  new_sample_values[i][Y] - old_sample_values[i][Y];
+        }
+        Vector column(m_model.size());
+        Matrix delta(m_model.size(), 2, 0.0);
+        VectorView deltax = delta.column_view(X);
+        VectorView deltay = delta.column_view(Y);
+        for (size_t i = 0; i < total_samples(); ++i)
+        {
+            if (diff(i, X) != 0)
+            {
+                column = m_psdinv_matrix->column_view(i);
+                column.scale(diff(i, X));
+                deltax += column;
+            }
+            if (diff(i, Y) != 0)
+            {
+                column = m_psdinv_matrix->column_view(i);
+                column.scale(diff(i, Y));
+                deltay += column;
+            }
+        }
+        m_solution += delta;
+        return m_solution;
+    }
+
+    solution_type& result()
+    {
+        return m_solution;
+    }
+
+private:
+    solution_type m_solution;
+
+}; // end class lsf_solution<ModelT, Point>
+
+
+
+
+template< typename ModelT,
+          bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
+class lsf_with_fixed_terms
+{
+};
+
+
+// fitting tool for completely unknown models
+template< typename ModelT>
+class lsf_with_fixed_terms<ModelT, false>
+    : public lsf_solution<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef typename model_type::parameter_type         parameter_type;
+    typedef typename model_type::value_type             value_type;
+    typedef lsf_solution<model_type>                    base_type;
+    typedef typename base_type::solution_type           solution_type;
+
+    using base_type::total_samples;
+    using base_type::is_full;
+    using base_type::m_matrix;
+    using base_type::m_total_samples;
+    using base_type::m_model;
+
+  public:
+      lsf_with_fixed_terms<ModelT, false>( model_type const& _model,
+                                           size_t forecasted_samples )
+        : base_type(_model, forecasted_samples)
+    {
+    }
+
+    void append(parameter_type const& sample_parameter)
+    {
+        assert(!is_full());
+        VectorView row = m_matrix.row_view(total_samples());
+        m_model.feed(row, sample_parameter);
+        ++m_total_samples;
+    }
+
+    void append_copy(size_t sample_index)
+    {
+        assert(!is_full());
+        assert(sample_index < total_samples());
+        VectorView dest_row = m_matrix.row_view(total_samples());
+        VectorView source_row = m_matrix.row_view(sample_index);
+        dest_row = source_row;
+        ++m_total_samples;
+    }
+
+}; // end class lsf_with_fixed_terms<ModelT, false>
+
+
+// fitting tool for partially known models
+template< typename ModelT>
+class lsf_with_fixed_terms<ModelT, true>
+    : public lsf_solution<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef typename model_type::parameter_type         parameter_type;
+    typedef typename model_type::value_type             value_type;
+    typedef lsf_solution<model_type>                    base_type;
+    typedef typename base_type::solution_type           solution_type;
+
+    using base_type::total_samples;
+    using base_type::is_full;
+    using base_type::m_matrix;
+    using base_type::m_total_samples;
+    using base_type::m_model;
+
+  public:
+    lsf_with_fixed_terms<ModelT, true>( model_type const& _model,
+                                        size_t forecasted_samples )
+        : base_type(_model, forecasted_samples),
+          m_vector(forecasted_samples),
+          m_vector_view(NULL)
+    {
+    }
+    void append(parameter_type const& sample_parameter)
+    {
+        assert(!is_full());
+        VectorView row = m_matrix.row_view(total_samples());
+        m_model.feed(row, m_vector[total_samples()], sample_parameter);
+        ++m_total_samples;
+    }
+
+    void append_copy(size_t sample_index)
+    {
+        assert(!is_full());
+        assert(sample_index < total_samples());
+        VectorView dest_row = m_matrix.row_view(total_samples());
+        VectorView source_row = m_matrix.row_view(sample_index);
+        dest_row = source_row;
+        m_vector[total_samples()] = m_vector[sample_index];
+        ++m_total_samples;
+    }
+
+    void update()
+    {
+        base_type::update();
+        if (total_samples() == 0) return;
+        if (m_vector_view != NULL)
+        {
+            delete m_vector_view;
+        }
+        m_vector_view = new VectorView(m_vector, base_type::total_samples());
+        assert(m_vector_view != NULL);
+    }
+
+    
+    ~lsf_with_fixed_terms<model_type, true>() override
+    {
+        if (m_vector_view != NULL)
+        {
+            delete m_vector_view;
+        }
+    }
+
+  protected:
+    Vector          m_vector;
+    VectorView*     m_vector_view;
+
+}; // end class lsf_with_fixed_terms<ModelT, true>
+
+
+} // end namespace detail
+
+
+
+
+template< typename ModelT,
+          typename ValueType = typename ModelT::value_type,
+          bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
+class least_squeares_fitter
+{
+};
+
+
+template< typename ModelT, typename ValueType >
+class least_squeares_fitter<ModelT, ValueType, false>
+    : public detail::lsf_with_fixed_terms<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef detail::lsf_with_fixed_terms<model_type>    base_type;
+    typedef typename base_type::parameter_type          parameter_type;
+    typedef typename base_type::value_type              value_type;
+    typedef typename base_type::solution_type           solution_type;
+
+  public:
+    least_squeares_fitter<ModelT, ValueType, false>( model_type const& _model,
+                                                     size_t forecasted_samples )
+          : base_type(_model, forecasted_samples)
+  {
+  }
+}; // end class least_squeares_fitter<ModelT, ValueType, true>
+
+
+template< typename ModelT>
+class least_squeares_fitter<ModelT, double, true>
+    : public detail::lsf_with_fixed_terms<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef detail::lsf_with_fixed_terms<model_type>    base_type;
+    typedef typename base_type::parameter_type          parameter_type;
+    typedef typename base_type::value_type              value_type;
+    typedef typename base_type::solution_type           solution_type;
+
+    using base_type::m_vector_view;
+    //using base_type::result; // VSC legacy support
+    solution_type& result( std::vector<Point> const& old_sample_values,
+                           std::vector<Point> const& new_sample_values )
+    {
+        return base_type::result(old_sample_values, new_sample_values);
+    }
+ 
+    solution_type& result()
+    {
+        return base_type::result();
+    }
+
+  public:
+    least_squeares_fitter<ModelT, double, true>( model_type const& _model,
+                                                 size_t forecasted_samples )
+        : base_type(_model, forecasted_samples)
+    {
+    }
+
+    template< typename VectorT >
+    solution_type& result(VectorT const& sample_values)
+    {
+        assert(sample_values.size() == m_vector_view->size());
+        Vector sv(sample_values.size());
+        for (size_t i = 0; i < sv.size(); ++i)
+            sv[i] = sample_values[i] - (*m_vector_view)[i];
+        return base_type::result(sv);
+    }
+
+}; // end class least_squeares_fitter<ModelT, double, true>
+
+
+template< typename ModelT>
+class least_squeares_fitter<ModelT, Point, true>
+    : public detail::lsf_with_fixed_terms<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef detail::lsf_with_fixed_terms<model_type>    base_type;
+    typedef typename base_type::parameter_type          parameter_type;
+    typedef typename base_type::value_type              value_type;
+    typedef typename base_type::solution_type           solution_type;
+
+    using base_type::m_vector_view;
+    //using base_type::result; // VCS legacy support
+    solution_type& result( std::vector<Point> const& old_sample_values,
+                           std::vector<Point> const& new_sample_values )
+    {
+        return base_type::result(old_sample_values, new_sample_values);
+    }
+ 
+    solution_type& result()
+    {
+        return base_type::result();
+    }
+
+
+  public:
+    least_squeares_fitter<ModelT, Point, true>( model_type const& _model,
+                                                size_t forecasted_samples )
+        : base_type(_model, forecasted_samples)
+    {
+    }
+
+    solution_type& result(std::vector<Point> const& sample_values)
+    {
+        assert(sample_values.size() == m_vector_view->size());
+        NL::Matrix sv(sample_values.size(), 2);
+        for (size_t i = 0; i < sample_values.size(); ++i)
+        {
+            sv(i, X) = sample_values[i][X] - (*m_vector_view)[i];
+            sv(i, Y) = sample_values[i][Y] - (*m_vector_view)[i];
+        }
+        return base_type::result(sv);
+    }
+
+};  // end class least_squeares_fitter<ModelT, Point, true>
+
+
+}  // end namespace NL
+}  // end namespace Geom
+
+
+
+#endif  // _NL_FITTING_TOOL_H_
+
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :