update 2geom (svn rev1433)

(bzr r6134)

6 files changed, 2328 insertions, 481 deletions

diff --git a/src/2geom/numeric/fitting-model.h b/src/2geom/numeric/fitting-model.h
new file mode 100644
index 000000000..145be40e4
--- /dev/null
+++ b/src/2geom/numeric/fitting-model.h
@@ -0,0 +1,423 @@
+/*

+ * Fitting Models for Geom Types

+ *

+ * 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_MODEL_H_

+#define _NL_FITTING_MODEL_H_

+

+

+#include <2geom/d2.h>

+#include <2geom/sbasis.h>

+#include <2geom/bezier.h>

+#include <2geom/bezier-curve.h>

+#include <2geom/poly.h>

+#include <2geom/ellipse.h>

+#include <2geom/utils.h>

+

+

+namespace Geom { namespace NL {

+

+

+/*

+ *   completely unknown models must inherit from this template class;

+ *   example: the model a*x^2 + b*x + c = 0 to be solved wrt a, b, c;

+ *   example: the model A(t) = known_sample_value_at(t) to be solved wrt

+ *       the coefficients of the curve A(t) expressed in S-Basis form;

+ *   parameter type: the type of x and t variable in the examples above;

+ *   value type:     the type of the known sample values (in the first example

+ *                   is constant )

+ *   instance type:  the type of the objects produced by using

+ *                   the fitting raw data solution

+ */

+template< typename ParameterType, typename ValueType, typename InstanceType >

+class LinearFittingModel

+{

+  public:

+    typedef ParameterType       parameter_type;

+    typedef ValueType           value_type;

+    typedef InstanceType        instance_type;

+

+    static const bool WITH_FIXED_TERMS = false;

+

+    /*

+     * a LinearFittingModel must implement the following methods:

+     *

+     * void feed( VectorView & vector,

+     *            parameter_type const& sample_parameter ) const;

+     *

+     * size_t size() const;

+     *

+     * void instance(instance_type &, raw_type const& raw_data) const;

+     *

+     */

+};

+

+

+/*

+ *   partially known models must inherit from this template class

+ *   example: the model a*x^2 + 2*x + c = 0 to be solved wrt a and c

+ */

+template< typename ParameterType, typename ValueType, typename InstanceType >

+class LinearFittingModelWithFixedTerms

+{

+  public:

+    typedef ParameterType       parameter_type;

+    typedef ValueType           value_type;

+    typedef InstanceType        instance_type;

+

+    static const bool WITH_FIXED_TERMS = true;

+

+    /*

+     * a LinearFittingModelWithFixedTerms must implement the following methods:

+     *

+     * void feed( VectorView & vector,

+     *            value_type & fixed_term,

+     *            parameter_type const& sample_parameter ) const;

+     *

+     * size_t size() const;

+     *

+     * void instance(instance_type &, raw_type const& raw_data) const;

+     *

+     */

+

+

+};

+

+

+// incomplete model, it can be inherited to make up different kinds of

+// instance type; the raw data is a vector of coefficients of a polynomial

+// rapresented in standard power basis

+template< typename InstanceType >

+class LFMPowerBasis

+    : public LinearFittingModel<double, double, InstanceType>

+{

+  public:

+    LFMPowerBasis(size_t degree)

+        : m_size(degree + 1)

+    {

+    }

+

+    void feed( VectorView & coeff, double sample_parameter ) const

+    {

+        coeff[0] = 1;

+        double x_i = 1;

+        for (size_t i = 1; i < coeff.size(); ++i)

+        {

+          x_i *= sample_parameter;

+          coeff[i] = x_i;

+        }

+    }

+

+    size_t size() const

+    {

+        return m_size;

+    }

+

+  private:

+    size_t m_size;

+};

+

+

+// this model generates Geom::Poly objects

+class LFMPoly

+    : public LFMPowerBasis<Poly>

+{

+  public:

+    LFMPoly(size_t degree)

+        : LFMPowerBasis<Poly>(degree)

+    {

+    }

+

+    void instance(Poly & poly, ConstVectorView const& raw_data) const

+    {

+        poly.clear();

+        poly.resize(size());

+        for (size_t i = 0; i < raw_data.size(); ++i)

+        {

+            poly[i] =  raw_data[i];

+        }

+    }

+};

+

+

+// incomplete model, it can be inherited to make up different kinds of

+// instance type; the raw data is a vector of coefficients of a polynomial

+// rapresented in standard power basis with leading term coefficient equal to 1

+template< typename InstanceType >

+class LFMNormalizedPowerBasis

+    : public LinearFittingModelWithFixedTerms<double, double, InstanceType>

+{

+  public:

+    LFMNormalizedPowerBasis(size_t _degree)

+        : m_model( _degree - 1)

+    {

+        assert(_degree > 0);

+    }

+

+

+    void feed( VectorView & coeff,

+               double & known_term,

+               double sample_parameter ) const

+    {

+        m_model.feed(coeff, sample_parameter);

+        known_term = coeff[m_model.size()-1] * sample_parameter;

+    }

+

+    size_t size() const

+    {

+        return m_model.size();

+    }

+

+  private:

+    LFMPowerBasis<InstanceType> m_model;

+};

+

+

+// incomplete model, it can be inherited to make up different kinds of

+// instance type; the raw data is a vector of coefficients of the equation

+// of an ellipse curve

+template< typename InstanceType >

+class LFMEllipseEquation

+    : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>

+{

+  public:

+    void feed( VectorView & coeff, double & fixed_term, Point const& p ) const

+    {

+        coeff[0] = p[X] * p[Y];

+        coeff[1] = p[Y] * p[Y];

+        coeff[2] = p[X];

+        coeff[3] = p[Y];

+        coeff[4] = 1;

+        fixed_term = p[X] * p[X];

+    }

+

+    size_t size() const

+    {

+        return 5;

+    }

+};

+

+

+// this model generates Ellipse curves

+class LFMEllipse

+    : public LFMEllipseEquation<Ellipse>

+{

+  public:

+    void instance(Ellipse & e, ConstVectorView const& coeff) const

+    {

+        e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);

+    }

+};

+

+

+// this model generates SBasis objects

+class LFMSBasis

+    : public LinearFittingModel<double, double, SBasis>

+{

+  public:

+    LFMSBasis( size_t _order )

+        : m_size( 2*(_order+1) ),

+          m_order(_order)

+    {

+    }

+

+    void feed( VectorView & coeff, double t ) const

+    {

+        double u0 = 1-t;

+        double u1 = t;

+        double s = u0 * u1;

+        coeff[0] = u0;

+        coeff[1] = u1;

+        for (size_t i = 2; i < size(); i+=2)

+        {

+            u0 *= s;

+            u1 *= s;

+            coeff[i] = u0;

+            coeff[i+1] = u1;

+        }

+    }

+

+    size_t size() const

+    {

+        return m_size;

+    }

+

+    void instance(SBasis & sb, ConstVectorView const& raw_data) const

+    {

+        sb.clear();

+        sb.resize(m_order+1);

+        for (unsigned int i = 0, k = 0; i < raw_data.size(); i+=2, ++k)

+        {

+            sb[k][0] = raw_data[i];

+            sb[k][1] = raw_data[i+1];

+        }

+    }

+

+  private:

+    size_t m_size;

+    size_t m_order;

+};

+

+

+// this model generates D2<SBasis> objects

+class LFMD2SBasis

+    : public LinearFittingModel< double, Point, D2<SBasis> >

+{

+  public:

+    LFMD2SBasis( size_t _order )

+        : mosb(_order)

+    {

+    }

+

+    void feed( VectorView & coeff, double t ) const

+    {

+        mosb.feed(coeff, t);

+    }

+

+    size_t size() const

+    {

+        return mosb.size();

+    }

+

+    void instance(D2<SBasis> & d2sb, ConstMatrixView const& raw_data) const

+    {

+        mosb.instance(d2sb[X], raw_data.column_const_view(X));

+        mosb.instance(d2sb[Y], raw_data.column_const_view(Y));

+    }

+

+  private:

+    LFMSBasis mosb;

+};

+

+

+// this model generates Bezier objects

+class LFMBezier

+    : public LinearFittingModel<double, double, Bezier>

+{

+  public:

+    LFMBezier( size_t _order )

+        : m_size(_order + 1),

+          m_order(_order)

+    {

+        binomial_coefficients(m_bc, m_order);

+    }

+

+    void feed( VectorView & coeff, double t ) const

+    {

+        double s = 1;

+        for (size_t i = 0; i < size(); ++i)

+        {

+            coeff[i] = s * m_bc[i];

+            s *= t;

+        }

+        double u = 1-t;

+        s = 1;

+        for (size_t i = size()-1; i > 0; --i)

+        {

+            coeff[i] *= s;

+            s *= u;

+        }

+        coeff[0] *= s;

+    }

+

+    size_t size() const

+    {

+        return m_size;

+    }

+

+    void instance(Bezier & b, ConstVectorView const& raw_data) const

+    {

+        assert(b.size() == raw_data.size());

+        for (unsigned int i = 0; i < raw_data.size(); ++i)

+        {

+            b[i] = raw_data[i];

+        }

+    }

+

+  private:

+    size_t m_size;

+    size_t m_order;

+    std::vector<size_t> m_bc;

+};

+

+

+// this model generates Bezier curves

+template< unsigned int N >

+class LFMBezierCurve

+    : public LinearFittingModel< double, Point, BezierCurve<N> >

+{

+  public:

+    LFMBezierCurve( size_t _order )

+        : mob(_order)

+    {

+    }

+

+    void feed( VectorView & coeff, double t ) const

+    {

+        mob.feed(coeff, t);

+    }

+

+    size_t size() const

+    {

+        return mob.size();

+    }

+

+    void instance(BezierCurve<N> & bc, ConstMatrixView const& raw_data) const

+    {

+        Bezier bx(size()-1);

+        Bezier by(size()-1);

+        mob.instance(bx, raw_data.column_const_view(X));

+        mob.instance(by, raw_data.column_const_view(Y));

+        bc = BezierCurve<N>(bx, by);

+    }

+

+  private:

+    LFMBezier mob;

+};

+

+}  // end namespace NL

+}  // end namespace Geom

+

+

+#endif // _NL_FITTING_MODEL_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:encoding=utf-8:textwidth=99 :

diff --git a/src/2geom/numeric/fitting-tool.h b/src/2geom/numeric/fitting-tool.h
new file mode 100644
index 000000000..edacc663a
--- /dev/null
+++ b/src/2geom/numeric/fitting-tool.h
@@ -0,0 +1,532 @@
+/*

+ * 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>

+

+

+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);

+    }

+

+    virtual

+    ~lsf_with_fixed_terms<model_type, true>()

+    {

+        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;

+

+  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;

+

+  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:encoding=utf-8:textwidth=99 :

diff --git a/src/2geom/numeric/linear_system.h b/src/2geom/numeric/linear_system.h
index 7699c5224..5b516c9e6 100644
--- a/src/2geom/numeric/linear_system.h
+++ b/src/2geom/numeric/linear_system.h
@@ -1,89 +1,138 @@
-#ifndef _NL_LINEAR_SYSTEM_H_
-#define _NL_LINEAR_SYSTEM_H_
-
-
-#include <cassert>
-
-#include <gsl/gsl_linalg.h>
-
-#include "2geom/numeric/matrix.h"
-#include "2geom/numeric/vector.h"
-
-
-namespace Geom { namespace NL {
-
-
-class LinearSystem
-{
-public:
-	LinearSystem(Matrix & _matrix, Vector & _vector)
-		: m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())
-	{
-	}
-	
-	const Vector & LU_solve()
-	{
-		assert( matrix().rows() == matrix().columns() 
-				&& matrix().rows() == vector().size() );
-		int s;
-		gsl_permutation * p = gsl_permutation_alloc(matrix().rows());
-		gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s);
-		gsl_linalg_LU_solve( matrix().get_gsl_matrix(), 
-							 p, 
-				             vector().get_gsl_vector(), 
-				             m_solution.get_gsl_vector()
-				           );
-		gsl_permutation_free(p);
-		return solution();
-	}
-	
-	const Vector & SV_solve()
-	{
-		assert( matrix().rows() >= matrix().columns()
-				&& matrix().rows() == vector().size() );
-		
-		gsl_matrix* U = matrix().get_gsl_matrix();
-		gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns());
-		gsl_vector* S = gsl_vector_alloc(matrix().columns());
-		gsl_vector* work = gsl_vector_alloc(matrix().columns());
-		
-		gsl_linalg_SV_decomp( U, V, S, work );
-		
-		gsl_vector* b = vector().get_gsl_vector();
-		gsl_vector* x = m_solution.get_gsl_vector();
-		
-		gsl_linalg_SV_solve( U, V, S, b, x);
-		
-		gsl_matrix_free(V);
-		gsl_vector_free(S);
-		gsl_vector_free(work);
-		
-		return solution();			  
-	}
-	
-	Matrix & matrix()
-	{
-		return m_matrix;
-	}
-	
-	Vector & vector()
-	{
-		return m_vector;
-	}
-	
-	const Vector & solution() const
-	{
-		return m_solution;
-	}
-	
-private:
-	Matrix & m_matrix;
-	Vector & m_vector;
-	Vector m_solution;
-};
-
-
-} } // end namespaces
-
-
-#endif /*_NL_LINEAR_SYSTEM_H_*/
+/*

+ * LinearSystem class wraps some gsl routines for solving linear systems

+ *

+ * 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_LINEAR_SYSTEM_H_

+#define _NL_LINEAR_SYSTEM_H_

+

+

+#include <cassert>

+

+#include <gsl/gsl_linalg.h>

+

+#include <2geom/numeric/matrix.h>

+#include <2geom/numeric/vector.h>

+

+

+namespace Geom { namespace NL {

+

+

+class LinearSystem

+{

+public:

+	LinearSystem(MatrixView & _matrix, VectorView & _vector)

+		: m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())

+	{

+	}

+	

+	LinearSystem(Matrix & _matrix, Vector & _vector)

+		: m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())

+	{

+	}

+	

+	const Vector & LU_solve()

+	{

+		assert( matrix().rows() == matrix().columns() 

+				&& matrix().rows() == vector().size() );

+		int s;

+		gsl_permutation * p = gsl_permutation_alloc(matrix().rows());

+		gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s);

+		gsl_linalg_LU_solve( matrix().get_gsl_matrix(), 

+							 p, 

+				             vector().get_gsl_vector(), 

+				             m_solution.get_gsl_vector()

+				           );

+		gsl_permutation_free(p);

+		return solution();

+	}

+	

+	const Vector & SV_solve()

+	{

+		assert( matrix().rows() >= matrix().columns()

+				&& matrix().rows() == vector().size() );

+		

+		gsl_matrix* U = matrix().get_gsl_matrix();

+		gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns());

+		gsl_vector* S = gsl_vector_alloc(matrix().columns());

+		gsl_vector* work = gsl_vector_alloc(matrix().columns());

+		

+		gsl_linalg_SV_decomp( U, V, S, work );

+		

+		gsl_vector* b = vector().get_gsl_vector();

+		gsl_vector* x = m_solution.get_gsl_vector();

+		

+		gsl_linalg_SV_solve( U, V, S, b, x);

+		

+		gsl_matrix_free(V);

+		gsl_vector_free(S);

+		gsl_vector_free(work);

+		

+		return solution();			  

+	}

+	

+	MatrixView & matrix()

+	{

+		return m_matrix;

+	}

+	

+	VectorView & vector()

+	{

+		return m_vector;

+	}

+	

+	const Vector & solution() const

+	{

+		return m_solution;

+	}

+	

+private:

+	MatrixView m_matrix;

+	VectorView m_vector;

+	Vector m_solution;

+};

+

+

+} } // end namespaces

+

+

+#endif /*_NL_LINEAR_SYSTEM_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:encoding=utf-8:textwidth=99 :

diff --git a/src/2geom/numeric/matrix.cpp b/src/2geom/numeric/matrix.cpp
new file mode 100644
index 000000000..dde80dfb6
--- /dev/null
+++ b/src/2geom/numeric/matrix.cpp
@@ -0,0 +1,115 @@
+/*
+ * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines;
+ * "views" mimic the semantic of C++ references: any operation performed
+ * on a "view" is actually performed on the "viewed object"
+ *
+ * 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.
+ */
+
+
+#include <numeric/matrix.h>
+#include <numeric/vector.h>
+
+
+
+namespace Geom { namespace NL {
+
+
+Vector operator*( detail::BaseMatrixImpl const& A,
+                  detail::BaseVectorImpl const& v )
+{
+    assert(A.columns() == v.size());
+
+    Vector result(A.rows(), 0.0);
+    for (size_t i = 0; i < A.rows(); ++i)
+        for (size_t j = 0; j < A.columns(); ++j)
+            result[i] += A(i,j) * v[j];
+
+    return result;
+}
+
+Matrix operator*( detail::BaseMatrixImpl const& A,
+                  detail::BaseMatrixImpl const& B )
+{
+    assert(A.columns() == B.rows());
+
+    Matrix C(A.rows(), B.columns(), 0.0);
+    for (size_t i = 0; i < C.rows(); ++i)
+        for (size_t j = 0; j < C.columns(); ++j)
+            for (size_t k = 0; k < A.columns(); ++k)
+                C(i,j) += A(i,k) * B(k, j);
+
+    return C;
+}
+
+Matrix pseudo_inverse(detail::BaseMatrixImpl const& A)
+{
+
+    Matrix U(A);
+    Matrix V(A.columns(), A.columns());
+    Vector s(A.columns());
+    gsl_vector* work = gsl_vector_alloc(A.columns());
+
+    gsl_linalg_SV_decomp( U.get_gsl_matrix(),
+                          V.get_gsl_matrix(),
+                          s.get_gsl_vector(),
+                          work );
+
+    Matrix P(A.columns(), A.rows(), 0.0);
+
+    int sz = s.size();
+    while ( sz-- > 0 && s[sz] == 0 ) {}
+    ++sz;
+    if (sz == 0) return P;
+    VectorView sv(s, sz);
+
+    for (size_t i = 0; i < sv.size(); ++i)
+    {
+        VectorView v = V.column_view(i);
+        v.scale(1/sv[i]);
+        for (size_t h = 0; h < P.rows(); ++h)
+            for (size_t k = 0; k < P.columns(); ++k)
+                P(h,k) += V(h,i) * U(k,i);
+    }
+
+    return P;
+}
+
+} }  // end namespaces
+
+/*
+  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:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/matrix.h b/src/2geom/numeric/matrix.h
index bdd647e53..64557a6f1 100644
--- a/src/2geom/numeric/matrix.h
+++ b/src/2geom/numeric/matrix.h
@@ -1,207 +1,555 @@
-#ifndef _NL_MATRIX_H_
-#define _NL_MATRIX_H_
-
-
-#include <cassert>
-#include <utility>
-#include <algorithm>
-
-#include <gsl/gsl_matrix.h>
-
-
-
-namespace Geom { namespace NL {
-
-class Matrix;
-void swap(Matrix & m1, Matrix & m2);
-
-
-class Matrix
-{
-public:
-	// the matrix is not inizialized
-	Matrix(size_t n1, size_t n2)
-		: m_rows(n1), m_columns(n2)
-	{
-		m_matrix = gsl_matrix_alloc(n1, n2);
-	}
-	
-	Matrix(size_t n1, size_t n2, double x)
-		:m_rows(n1), m_columns(n2)
-	{
-		m_matrix = gsl_matrix_alloc(n1, n2);
-		gsl_matrix_set_all(m_matrix, x);
-	}
-	
-	Matrix( Matrix const& _matrix )
-		: m_rows(_matrix.rows()), m_columns(_matrix.columns())
-	{
-		m_matrix = gsl_matrix_alloc(rows(), columns());
-		gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);
-	}
-	
-//	explicit
-//	Matrix( gsl_matrix* m, size_t n1, size_t n2)
-//		: m_rows(n1), m_columns(n2), m_matrix(m)
-//	{	
-//	}
-	
-	Matrix & operator=(Matrix const& _matrix)
-	{
-		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );
-		gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);
-		
-		return *this;
-	}
-	
-	virtual ~Matrix()
-	{
-		gsl_matrix_free(m_matrix);
-	}
-	
-	void set_all( double x = 0)
-	{
-		gsl_matrix_set_all(m_matrix, x);
-	}
-	
-	void set_identity()
-	{
-		gsl_matrix_set_identity(m_matrix);
-	}
-	
-	const double & operator() (size_t i, size_t j) const
-	{
-		return *gsl_matrix_const_ptr(m_matrix, i, j);
-	}
-	
-	double & operator() (size_t i, size_t j)
-	{
-		return *gsl_matrix_ptr(m_matrix, i, j);
-	}
-	
-	gsl_matrix* get_gsl_matrix()
-	{
-		return m_matrix;
-	}
-	
-	void swap_rows(size_t i, size_t j)
-	{
-		 gsl_matrix_swap_rows(m_matrix, i, j);
-	}
-	
-	void swap_columns(size_t i, size_t j)
-	{
-		gsl_matrix_swap_columns(m_matrix, i, j);
-	}
-	
-	Matrix & transpose()
-	{
-		gsl_matrix_transpose(m_matrix);
-		return (*this);
-	}
-	
-	Matrix & scale(double x)
-	{
-		gsl_matrix_scale(m_matrix, x);
-		return (*this);
-	}
-	
-	Matrix & translate(double x)
-	{
-		 gsl_matrix_add_constant(m_matrix, x);
-		 return (*this);
-	}
-	
-	Matrix & operator+=(Matrix const& _matrix)
-	{
-		gsl_matrix_add(m_matrix, _matrix.m_matrix);
-		return (*this);
-	}
-	
-	Matrix & operator-=(Matrix const& _matrix)
-	{
-		gsl_matrix_sub(m_matrix, _matrix.m_matrix);
-		return (*this);
-	}
-	
-	bool is_zero() const
-	{
-		return gsl_matrix_isnull(m_matrix);
-	}
-	
-	bool is_positive() const
-	{
-		return gsl_matrix_ispos(m_matrix);
-	}
-	
-	bool is_negative() const
-	{
-		return gsl_matrix_isneg(m_matrix);
-	}
-	
-	bool is_non_negative() const
-	{
-		for ( unsigned int i = 0; i < rows(); ++i )
-		{
-			for ( unsigned int j = 0; j < columns(); ++j )
-			{
-				if ( (*this)(i,j) < 0 ) return false;
-			}
-		}
-		return true;
-	}
-	
-	double max() const
-	{
-		return gsl_matrix_max(m_matrix);
-	}
-	
-	double min() const
-	{
-		return gsl_matrix_min(m_matrix);
-	}
-	
-	std::pair<size_t, size_t>
-	max_index() const
-	{
-		std::pair<size_t, size_t> indices;
-		gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));
-		return indices;
-	}
-	
-	std::pair<size_t, size_t>
-	min_index() const
-	{
-		std::pair<size_t, size_t> indices;
-		gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));
-		return indices;
-	}
-	
-	friend 
-	void swap(Matrix & m1, Matrix & m2);
-	
-	size_t rows() const
-	{
-		return m_rows;
-	}
-	
-	size_t columns() const
-	{
-		return m_columns;
-	}
-	
-private:
-	size_t m_rows, m_columns;
-	gsl_matrix* m_matrix;
-};
-
-void swap(Matrix & m1, Matrix & m2)
-{
-	assert( m1.rows() == m2.rows() && m1.columns() ==  m2.columns() );
-	std::swap(m1.m_matrix, m2.m_matrix);
-}
-
-
-} } // end namespaces
-
-#endif /*_NL_MATRIX_H_*/
+/*

+ * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines;

+ * "views" mimic the semantic of C++ references: any operation performed

+ * on a "view" is actually performed on the "viewed object"

+ *

+ * 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_MATRIX_H_

+#define _NL_MATRIX_H_

+

+#include <2geom/numeric/vector.h>

+

+#include <cassert>

+#include <utility>    // for std::pair

+#include <algorithm>  // for std::swap

+#include <sstream>

+#include <string>

+

+#include <gsl/gsl_matrix.h>

+#include <gsl/gsl_linalg.h>

+

+

+namespace Geom { namespace NL {

+

+namespace detail

+{

+

+class BaseMatrixImpl

+{

+  public:

+	virtual ~BaseMatrixImpl()

+	{

+	}

+

+	ConstVectorView row_const_view(size_t i) const

+	{

+		return ConstVectorView(gsl_matrix_const_row(m_matrix, i));

+	}

+

+	ConstVectorView column_const_view(size_t i) const

+	{

+		return ConstVectorView(gsl_matrix_const_column(m_matrix, i));

+	}

+

+	const double & operator() (size_t i, size_t j) const

+	{

+		return *gsl_matrix_const_ptr(m_matrix, i, j);

+	}

+

+	const gsl_matrix* get_gsl_matrix() const

+	{

+		return m_matrix;

+	}

+

+	bool is_zero() const

+	{

+		return gsl_matrix_isnull(m_matrix);

+	}

+

+	bool is_positive() const

+	{

+		return gsl_matrix_ispos(m_matrix);

+	}

+

+	bool is_negative() const

+	{

+		return gsl_matrix_isneg(m_matrix);

+	}

+

+	bool is_non_negative() const

+	{

+		for ( unsigned int i = 0; i < rows(); ++i )

+		{

+			for ( unsigned int j = 0; j < columns(); ++j )

+			{

+				if ( (*this)(i,j) < 0 ) return false;

+			}

+		}

+		return true;

+	}

+

+	double max() const

+	{

+		return gsl_matrix_max(m_matrix);

+	}

+

+	double min() const

+	{

+		return gsl_matrix_min(m_matrix);

+	}

+

+	std::pair<size_t, size_t>

+	max_index() const

+	{

+		std::pair<size_t, size_t> indices;

+		gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));

+		return indices;

+	}

+

+	std::pair<size_t, size_t>

+	min_index() const

+	{

+		std::pair<size_t, size_t> indices;

+		gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));

+		return indices;

+	}

+

+	size_t rows() const

+	{

+		return m_rows;

+	}

+

+	size_t columns() const

+	{

+		return m_columns;

+	}

+

+	std::string str() const;

+

+  protected:

+	size_t m_rows, m_columns;

+	gsl_matrix* m_matrix;

+

+};  // end class BaseMatrixImpl

+

+

+inline

+bool operator== (BaseMatrixImpl const& m1, BaseMatrixImpl const& m2)

+{

+	if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) return false;

+

+	for (size_t i = 0; i < m1.rows(); ++i)

+		for (size_t j = 0; j < m1.columns(); ++j)

+			if (m1(i,j) != m2(i,j)) return false;

+

+	return true;

+}

+

+template< class charT >

+inline

+std::basic_ostream<charT> &

+operator<< (std::basic_ostream<charT> & os, const BaseMatrixImpl & _matrix)

+{

+	if (_matrix.rows() == 0 || _matrix.columns() == 0) return os;

+

+	os << "[[" << _matrix(0,0);

+	for (size_t j = 1; j < _matrix.columns(); ++j)

+	{

+		os << ", " << _matrix(0,j);

+	}

+	os << "]";

+

+	for (size_t i = 1; i < _matrix.rows(); ++i)

+	{

+		os << ", [" << _matrix(i,0);

+		for (size_t j = 1; j < _matrix.columns(); ++j)

+		{

+			os << ", " << _matrix(i,j);

+		}

+		os << "]";

+	}

+	os << "]";

+	return os;

+}

+

+inline

+std::string BaseMatrixImpl::str() const

+{

+	std::ostringstream oss;

+	oss << (*this);

+	return oss.str();

+}

+

+

+class MatrixImpl : public BaseMatrixImpl

+{

+  public:

+

+	typedef BaseMatrixImpl base_type;

+

+	void set_all( double x )

+	{

+		gsl_matrix_set_all(m_matrix, x);

+	}

+

+	void set_identity()

+	{

+		gsl_matrix_set_identity(m_matrix);

+	}

+

+	using base_type::operator();

+

+	double & operator() (size_t i, size_t j)

+	{

+		return *gsl_matrix_ptr(m_matrix, i, j);

+	}

+

+	using base_type::get_gsl_matrix;

+

+	gsl_matrix* get_gsl_matrix()

+	{

+		return m_matrix;

+	}

+

+	VectorView row_view(size_t i)

+	{

+		return VectorView(gsl_matrix_row(m_matrix, i));

+	}

+

+	VectorView column_view(size_t i)

+	{

+		return VectorView(gsl_matrix_column(m_matrix, i));

+	}

+

+	void swap_rows(size_t i, size_t j)

+	{

+		 gsl_matrix_swap_rows(m_matrix, i, j);

+	}

+

+	void swap_columns(size_t i, size_t j)

+	{

+		gsl_matrix_swap_columns(m_matrix, i, j);

+	}

+

+	MatrixImpl & transpose()

+	{

+	    assert(columns() == rows());

+		gsl_matrix_transpose(m_matrix);

+		return (*this);

+	}

+

+	MatrixImpl & scale(double x)

+	{

+		gsl_matrix_scale(m_matrix, x);

+		return (*this);

+	}

+

+	MatrixImpl & translate(double x)

+	{

+		 gsl_matrix_add_constant(m_matrix, x);

+		 return (*this);

+	}

+

+	MatrixImpl & operator+=(base_type const& _matrix)

+	{

+		gsl_matrix_add(m_matrix, _matrix.get_gsl_matrix());

+		return (*this);

+	}

+

+	MatrixImpl & operator-=(base_type const& _matrix)

+	{

+		gsl_matrix_sub(m_matrix, _matrix.get_gsl_matrix());

+		return (*this);

+	}

+

+}; // end class MatrixImpl

+

+}  // end namespace detail

+

+

+using detail::operator==;

+using detail::operator<<;

+

+

+

+

+class Matrix: public detail::MatrixImpl

+{

+  public:

+	typedef detail::MatrixImpl base_type;

+

+  public:

+	// the matrix is not inizialized

+	Matrix(size_t n1, size_t n2)

+	{

+		m_rows = n1;

+		m_columns = n2;

+		m_matrix = gsl_matrix_alloc(n1, n2);

+	}

+

+	Matrix(size_t n1, size_t n2, double x)

+	{

+		m_rows = n1;

+		m_columns = n2;

+		m_matrix = gsl_matrix_alloc(n1, n2);

+		gsl_matrix_set_all(m_matrix, x);

+	}

+

+	Matrix(Matrix const& _matrix)

+        : base_type()

+	{

+		m_rows = _matrix.rows();

+		m_columns = _matrix.columns();

+		m_matrix = gsl_matrix_alloc(rows(), columns());

+		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());

+	}

+

+	explicit

+	Matrix(base_type::base_type const& _matrix)

+	{

+		m_rows = _matrix.rows();

+		m_columns = _matrix.columns();

+		m_matrix = gsl_matrix_alloc(rows(), columns());

+		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());

+	}

+

+	Matrix & operator=(Matrix const& _matrix)

+	{

+		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );

+		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());

+		return *this;

+	}

+

+	Matrix & operator=(base_type::base_type const& _matrix)

+	{

+		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );

+		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());

+		return *this;

+	}

+

+	virtual ~Matrix()

+	{

+		gsl_matrix_free(m_matrix);

+	}

+

+	Matrix & transpose()

+	{

+		return static_cast<Matrix &>( base_type::transpose() );

+	}

+

+	Matrix & scale(double x)

+	{

+		return static_cast<Matrix &>( base_type::scale(x) );

+	}

+

+	Matrix & translate(double x)

+	{

+		return static_cast<Matrix &>( base_type::translate(x) );

+	}

+

+	Matrix & operator+=(base_type::base_type const& _matrix)

+	{

+		return static_cast<Matrix &>( base_type::operator+=(_matrix) );

+	}

+

+	Matrix & operator-=(base_type::base_type const& _matrix)

+	{

+		return static_cast<Matrix &>( base_type::operator-=(_matrix) );

+	}

+

+	friend

+	void swap(Matrix & m1, Matrix & m2);

+	friend

+	void swap_any(Matrix & m1, Matrix & m2);

+

+};  // end class Matrix

+

+

+// warning! this operation invalidates any view of the passed matrix objects

+inline

+void swap(Matrix & m1, Matrix & m2)

+{

+	assert( m1.rows() == m2.rows() && m1.columns() ==  m2.columns() );

+	std::swap(m1.m_matrix, m2.m_matrix);

+}

+

+inline

+void swap_any(Matrix & m1, Matrix & m2)

+{

+    std::swap(m1.m_matrix, m2.m_matrix);

+    std::swap(m1.m_rows, m2.m_rows);

+    std::swap(m1.m_columns, m2.m_columns);

+}

+

+

+

+class ConstMatrixView : public detail::BaseMatrixImpl

+{

+  public:

+	typedef detail::BaseMatrixImpl base_type;

+

+  public:

+	ConstMatrixView(const base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)

+		: m_matrix_view( gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2) )

+	{

+		m_rows = n1;

+		m_columns = n2;

+		m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );

+	}

+

+	ConstMatrixView(const ConstMatrixView & _matrix)

+		: base_type(),

+		  m_matrix_view(_matrix.m_matrix_view)

+	{

+		m_rows = _matrix.rows();

+		m_columns = _matrix.columns();

+		m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );

+	}

+

+	ConstMatrixView(const base_type & _matrix)

+		: m_matrix_view(gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), 0, 0, _matrix.rows(), _matrix.columns()))

+	{

+		m_rows = _matrix.rows();

+		m_columns = _matrix.columns();

+		m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );

+	}

+

+  private:

+	gsl_matrix_const_view m_matrix_view;

+

+};  // end class ConstMatrixView

+

+

+

+

+class MatrixView : public detail::MatrixImpl

+{

+  public:

+	typedef detail::MatrixImpl base_type;

+

+  public:

+	MatrixView(base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)

+	{

+		m_rows = n1;

+		m_columns = n2;

+		m_matrix_view

+			= gsl_matrix_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2);

+		m_matrix = &(m_matrix_view.matrix);

+	}

+

+	MatrixView(const MatrixView & _matrix)

+        : base_type()

+	{

+		m_rows = _matrix.rows();

+		m_columns = _matrix.columns();

+		m_matrix_view = _matrix.m_matrix_view;

+		m_matrix = &(m_matrix_view.matrix);

+	}

+

+	MatrixView(Matrix & _matrix)

+	{

+		m_rows = _matrix.rows();

+		m_columns = _matrix.columns();

+		m_matrix_view

+			= gsl_matrix_submatrix(_matrix.get_gsl_matrix(), 0, 0, rows(), columns());

+		m_matrix = &(m_matrix_view.matrix);

+	}

+

+	MatrixView & operator=(MatrixView const& _matrix)

+	{

+		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );

+		gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);

+		return *this;

+	}

+

+	MatrixView & operator=(base_type::base_type const& _matrix)

+	{

+		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );

+		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());

+		return *this;

+	}

+

+	MatrixView & transpose()

+	{

+		return static_cast<MatrixView &>( base_type::transpose() );

+	}

+

+	MatrixView & scale(double x)

+	{

+		return static_cast<MatrixView &>( base_type::scale(x) );

+	}

+

+	MatrixView & translate(double x)

+	{

+		return static_cast<MatrixView &>( base_type::translate(x) );

+	}

+

+	MatrixView & operator+=(base_type::base_type const& _matrix)

+	{

+		return static_cast<MatrixView &>( base_type::operator+=(_matrix) );

+	}

+

+	MatrixView & operator-=(base_type::base_type const& _matrix)

+	{

+		return static_cast<MatrixView &>( base_type::operator-=(_matrix) );

+	}

+

+	friend

+	void swap_view(MatrixView & m1, MatrixView & m2);

+

+  private:

+	gsl_matrix_view m_matrix_view;

+

+};  // end class MatrixView

+

+

+inline

+void swap_view(MatrixView & m1, MatrixView & m2)

+{

+	assert( m1.rows() == m2.rows() && m1.columns() ==  m2.columns() );

+	std::swap(m1.m_matrix_view, m2.m_matrix_view);

+}

+

+Vector operator*( detail::BaseMatrixImpl const& A,

+                  detail::BaseVectorImpl const& v );

+

+Matrix operator*( detail::BaseMatrixImpl const& A,

+                  detail::BaseMatrixImpl const& B );

+

+Matrix pseudo_inverse(detail::BaseMatrixImpl const& A);

+

+} } // end namespaces

+

+#endif /*_NL_MATRIX_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:encoding=utf-8:textwidth=99 :

diff --git a/src/2geom/numeric/vector.h b/src/2geom/numeric/vector.h
index b25861e22..43a39a1ac 100644
--- a/src/2geom/numeric/vector.h
+++ b/src/2geom/numeric/vector.h
@@ -1,185 +1,565 @@
-#ifndef _NL_VECTOR_H_
-#define _NL_VECTOR_H_
-
-#include <cassert>
-#include <utility>
-
-#include <gsl/gsl_vector.h>
-
-
-namespace Geom { namespace NL {
-
-class Vector;
-void swap(Vector & v1, Vector & v2);
-
-
-class Vector
-{
-public:
-	Vector(size_t n)
-		: m_size(n)
-	{
-		m_vector = gsl_vector_alloc(n);
-	}
-	
-	Vector(size_t n, double x)
-		: m_size(n)
-	{
-		m_vector = gsl_vector_alloc(n);
-		gsl_vector_set_all(m_vector, x);
-	}
-	
-	// create a vector with n elements all set to zero 
-	// but the i-th that is set to 1
-	Vector(size_t n, size_t i)
-		: m_size(n)
-	{
-		m_vector = gsl_vector_alloc(n);
-		gsl_vector_set_basis(m_vector, i);
-	}
-	
-	Vector(Vector const& _vector)
-		: m_size(_vector.size())
-	{
-		m_vector = gsl_vector_alloc(size());
-		gsl_vector_memcpy(m_vector, _vector.m_vector);
-	}
-	
-	virtual ~Vector()
-	{
-		gsl_vector_free(m_vector);
-	}
-	
-	Vector & operator=(Vector const& _vector)
-	{
-		assert( size() == _vector.size() );
-		gsl_vector_memcpy(m_vector, _vector.m_vector);
-		return (*this);
-	}
-	
-	void set_all(double x = 0)
-	{
-		gsl_vector_set_all(m_vector, x);
-	}
-
-	void set_basis(size_t i)
-	{
-		gsl_vector_set_basis(m_vector, i);
-	}
-	
-	double const& operator[](size_t i) const
-	{
-		return *gsl_vector_const_ptr(m_vector, i);
-	}
-	
-	double & operator[](size_t i)
-	{
-		return *gsl_vector_ptr(m_vector, i);
-	}
-	
-	gsl_vector* get_gsl_vector()
-	{
-		return m_vector;
-	}
-	
-	void swap_elements(size_t i, size_t j)
-	{
-		gsl_vector_swap_elements(m_vector, i, j);
-	}
-	
-	void reverse()
-	{
-		gsl_vector_reverse(m_vector);
-	}
-	
-	Vector & scale(double x)
-	{
-		gsl_vector_scale(m_vector, x);
-		return (*this);
-	}
-	
-	Vector & translate(double x)
-	{
-		gsl_vector_add_constant(m_vector, x);
-		return (*this);
-	}
-	
-	Vector & operator+=(Vector const& _vector)
-	{
-		gsl_vector_add(m_vector, _vector.m_vector);
-		return (*this);
-	}
-	
-	Vector & operator-=(Vector const& _vector)
-	{
-		gsl_vector_sub(m_vector, _vector.m_vector);
-		return (*this);
-	}
-	
-	bool is_zero() const
-	{
-		return gsl_vector_isnull(m_vector);
-	}
-	
-	bool is_positive() const
-	{
-		return gsl_vector_ispos(m_vector);
-	}
-	
-	bool is_negative() const
-	{
-		return gsl_vector_isneg(m_vector);
-	}
-	
-	bool is_non_negative() const
-	{
-		for ( size_t i = 0; i < size(); ++i )
-		{
-			if ( (*this)[i] < 0 ) return false;
-		}
-		return true;
-	}
-	
-	double max() const
-	{
-		return gsl_vector_max(m_vector);
-	}
-	
-	double min() const
-	{
-		return gsl_vector_min(m_vector);
-	}
-	
-	size_t max_index() const
-	{
-		return gsl_vector_max_index(m_vector);
-	}
-	
-	size_t min_index() const
-	{
-		return gsl_vector_min_index(m_vector);
-	}
-	
-	friend
-	void swap(Vector & v1, Vector & v2);
-	
-	size_t size() const
-	{
-		return m_size;
-	}
-	
-private:
-	size_t m_size;
-	gsl_vector* m_vector;
-};
-
-void swap(Vector & v1, Vector & v2)
-{
-	assert( v1.size() == v2.size() );
-	std::swap(v1.m_vector, v2.m_vector);
-}
-
-} } // end namespaces
-
-
-#endif /*_NL_VECTOR_H_*/
+/*

+ * Vector, VectorView, ConstVectorView classes wrap the gsl vector routines;

+ * "views" mimic the semantic of C++ references: any operation performed

+ * on a "view" is actually performed on the "viewed object"

+ *

+ * 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_VECTOR_H_

+#define _NL_VECTOR_H_

+

+#include <cassert>

+#include <algorithm> // for std::swap

+#include <vector>

+#include <sstream>

+#include <string>

+

+

+#include <gsl/gsl_vector.h>

+#include <gsl/gsl_blas.h>

+

+

+namespace Geom { namespace NL {

+

+namespace detail

+{

+

+class BaseVectorImpl

+{

+  public:

+	double const& operator[](size_t i) const

+	{

+		return *gsl_vector_const_ptr(m_vector, i);

+	}

+

+	const gsl_vector* get_gsl_vector() const

+	{

+		return m_vector;

+	}

+	bool is_zero() const

+	{

+		return gsl_vector_isnull(m_vector);

+	}

+

+	bool is_positive() const

+	{

+		return gsl_vector_ispos(m_vector);

+	}

+

+	bool is_negative() const

+	{

+		return gsl_vector_isneg(m_vector);

+	}

+

+	bool is_non_negative() const

+	{

+		for ( size_t i = 0; i < size(); ++i )

+		{

+			if ( (*this)[i] < 0 ) return false;

+		}

+		return true;

+	}

+

+	double max() const

+	{

+		return gsl_vector_max(m_vector);

+	}

+

+	double min() const

+	{

+		return gsl_vector_min(m_vector);

+	}

+

+	size_t max_index() const

+	{

+		return gsl_vector_max_index(m_vector);

+	}

+

+	size_t min_index() const

+	{

+		return gsl_vector_min_index(m_vector);

+	}

+

+	size_t size() const

+	{

+		return m_size;

+	}

+

+	std::string str() const;

+

+	virtual ~BaseVectorImpl()

+	{

+	}

+

+  protected:

+	size_t m_size;

+	gsl_vector* m_vector;

+

+};  // end class BaseVectorImpl

+

+

+inline

+bool operator== (BaseVectorImpl const& v1, BaseVectorImpl const& v2)

+{

+	if (v1.size() != v2.size())	return false;

+

+	for (size_t i = 0; i < v1.size(); ++i)

+	{

+		if (v1[i] != v2[i])  return false;

+	}

+	return true;

+}

+

+template< class charT >

+inline

+std::basic_ostream<charT> &

+operator<< (std::basic_ostream<charT> & os, const BaseVectorImpl & _vector)

+{

+	if (_vector.size() == 0 ) return os;

+	os << "[" << _vector[0];

+	for (unsigned int i = 1; i < _vector.size(); ++i)

+	{

+		os << ", " << _vector[i];

+	}

+	os << "]";

+	return os;

+}

+

+inline

+std::string BaseVectorImpl::str() const

+{

+	std::ostringstream oss;

+	oss << (*this);

+	return oss.str();

+}

+

+inline

+double dot(BaseVectorImpl const& v1, BaseVectorImpl const& v2)

+{

+    double result;

+    gsl_blas_ddot(v1.get_gsl_vector(), v2.get_gsl_vector(), &result);

+    return result;

+}

+

+

+class VectorImpl : public BaseVectorImpl

+{

+  public:

+	typedef BaseVectorImpl base_type;

+

+  public:

+	void set_all(double x)

+	{

+		gsl_vector_set_all(m_vector, x);

+	}

+

+	void set_basis(size_t i)

+	{

+		gsl_vector_set_basis(m_vector, i);

+	}

+

+	using base_type::operator[];

+

+	double & operator[](size_t i)

+	{

+		return *gsl_vector_ptr(m_vector, i);

+	}

+

+	using base_type::get_gsl_vector;

+

+	gsl_vector* get_gsl_vector()

+	{

+		return m_vector;

+	}

+

+	void swap_elements(size_t i, size_t j)

+	{

+		gsl_vector_swap_elements(m_vector, i, j);

+	}

+

+	void reverse()

+	{

+		gsl_vector_reverse(m_vector);

+	}

+

+	VectorImpl & scale(double x)

+	{

+		gsl_vector_scale(m_vector, x);

+		return (*this);

+	}

+

+	VectorImpl & translate(double x)

+	{

+		gsl_vector_add_constant(m_vector, x);

+		return (*this);

+	}

+

+	VectorImpl & operator+=(base_type const& _vector)

+	{

+		gsl_vector_add(m_vector, _vector.get_gsl_vector());

+		return (*this);

+	}

+

+	VectorImpl & operator-=(base_type const& _vector)

+	{

+		gsl_vector_sub(m_vector, _vector.get_gsl_vector());

+		return (*this);

+	}

+

+};  // end class VectorImpl

+

+}  // end namespace detail

+

+

+using detail::operator==;

+using detail::operator<<;

+

+class Vector : public detail::VectorImpl

+{

+  public:

+	typedef detail::VectorImpl base_type;

+

+  public:

+	Vector(size_t n)

+	{

+		m_size = n;

+		m_vector = gsl_vector_alloc(n);

+	}

+

+	Vector(size_t n, double x)

+	{

+		m_size = n;

+		m_vector = gsl_vector_alloc(n);

+		gsl_vector_set_all(m_vector, x);

+	}

+

+	// create a vector with n elements all set to zero

+	// but the i-th that is set to 1

+	Vector(size_t n, size_t i)

+	{

+		m_size = n;

+		m_vector = gsl_vector_alloc(n);

+		gsl_vector_set_basis(m_vector, i);

+	}

+

+	Vector(Vector const& _vector)

+        : base_type()

+	{

+		m_size = _vector.size();

+		m_vector = gsl_vector_alloc(size());

+		gsl_vector_memcpy(m_vector, _vector.m_vector);

+	}

+

+	explicit

+	Vector(base_type::base_type const& _vector)

+	{

+		m_size = _vector.size();

+		m_vector = gsl_vector_alloc(size());

+		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());

+	}

+

+	virtual ~Vector()

+	{

+		gsl_vector_free(m_vector);

+	}

+

+

+	Vector & operator=(Vector const& _vector)

+	{

+		assert( size() == _vector.size() );

+		gsl_vector_memcpy(m_vector, _vector.m_vector);

+		return (*this);

+	}

+

+	Vector & operator=(base_type::base_type const& _vector)

+	{

+		assert( size() == _vector.size() );

+		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());

+		return (*this);

+	}

+

+	Vector & scale(double x)

+	{

+		return static_cast<Vector&>( base_type::scale(x) );

+	}

+

+	Vector & translate(double x)

+	{

+		return static_cast<Vector&>( base_type::translate(x) );

+	}

+

+	Vector & operator+=(base_type::base_type const& _vector)

+	{

+		return static_cast<Vector&>( base_type::operator+=(_vector) );

+	}

+

+	Vector & operator-=(base_type::base_type const& _vector)

+	{

+		return static_cast<Vector&>( base_type::operator-=(_vector) );

+	}

+

+	friend

+	void swap(Vector & v1, Vector & v2);

+	friend

+	void swap_any(Vector & v1, Vector & v2);

+

+}; // end class Vector

+

+

+// warning! these operations invalidate any view of the passed vector objects

+inline

+void swap(Vector & v1, Vector & v2)

+{

+	assert( v1.size() == v2.size() );

+	std::swap(v1.m_vector, v2.m_vector);

+}

+

+inline

+void swap_any(Vector & v1, Vector & v2)

+{

+    std::swap(v1.m_vector, v2.m_vector);

+    std::swap(v1.m_size, v2.m_size);

+}

+

+

+class ConstVectorView : public detail::BaseVectorImpl

+{

+  public:

+	typedef detail::BaseVectorImpl base_type;

+

+  public:

+	ConstVectorView(const base_type & _vector, size_t n, size_t offset = 0)

+		: m_vector_view( gsl_vector_const_subvector(_vector.get_gsl_vector(), offset, n) )

+	{

+		m_size = n;

+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

+	}

+

+	ConstVectorView(const base_type & _vector, size_t n, size_t offset , size_t stride)

+		: m_vector_view( gsl_vector_const_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n) )

+	{

+		m_size = n;

+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

+	}

+

+    ConstVectorView(const double* _vector, size_t n, size_t offset = 0)

+        : m_vector_view( gsl_vector_const_view_array(_vector + offset, n) )

+    {

+        m_size = n;

+        m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

+    }

+

+    ConstVectorView(const double* _vector, size_t n, size_t offset, size_t stride)

+        : m_vector_view( gsl_vector_const_view_array_with_stride(_vector + offset, stride, n) )

+    {

+        m_size = n;

+        m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

+    }

+

+	explicit

+	ConstVectorView(gsl_vector_const_view  _gsl_vector_view)

+		: m_vector_view(_gsl_vector_view)

+	{

+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

+		m_size = m_vector->size;

+	}

+

+    explicit

+    ConstVectorView(const std::vector<double>&  _vector)

+        : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) )

+    {

+        m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

+        m_size = _vector.size();

+    }

+

+	ConstVectorView(const ConstVectorView & _vector)

+		: base_type(),

+		  m_vector_view(_vector.m_vector_view)

+	{

+		m_size = _vector.size();

+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

+	}

+

+	ConstVectorView(const base_type & _vector)

+		: m_vector_view(gsl_vector_const_subvector(_vector.get_gsl_vector(), 0, _vector.size()))

+	{

+		m_size = _vector.size();

+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

+	}

+

+  private:

+	gsl_vector_const_view m_vector_view;

+

+}; // end class  ConstVectorView

+

+

+

+

+class VectorView : public detail::VectorImpl

+{

+  public:

+	typedef detail::VectorImpl base_type;

+

+  public:

+	VectorView(base_type & _vector, size_t n, size_t offset = 0, size_t stride = 1)

+	{

+		m_size = n;

+		if (stride == 1)

+		{

+			m_vector_view

+				= gsl_vector_subvector(_vector.get_gsl_vector(), offset, n);

+			m_vector = &(m_vector_view.vector);

+		}

+		else

+		{

+			m_vector_view

+				= gsl_vector_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n);

+			m_vector = &(m_vector_view.vector);

+		}

+	}

+

+    VectorView(double* _vector, size_t n, size_t offset = 0, size_t stride = 1)

+    {

+        m_size = n;

+        if (stride == 1)

+        {

+            m_vector_view

+                = gsl_vector_view_array(_vector + offset, n);

+            m_vector = &(m_vector_view.vector);

+        }

+        else

+        {

+            m_vector_view

+                = gsl_vector_view_array_with_stride(_vector + offset, stride, n);

+            m_vector = &(m_vector_view.vector);

+        }

+

+    }

+

+	VectorView(const VectorView & _vector)

+        : base_type()

+	{

+		m_size = _vector.size();

+		m_vector_view = _vector.m_vector_view;

+		m_vector = &(m_vector_view.vector);

+	}

+

+	VectorView(Vector & _vector)

+	{

+		m_size = _vector.size();

+		m_vector_view = gsl_vector_subvector(_vector.get_gsl_vector(), 0, size());

+		m_vector = &(m_vector_view.vector);

+	}

+

+	explicit

+	VectorView(gsl_vector_view _gsl_vector_view)

+		: m_vector_view(_gsl_vector_view)

+	{

+		m_vector = &(m_vector_view.vector);

+		m_size = m_vector->size;

+	}

+

+	explicit

+	VectorView(std::vector<double> & _vector)

+	{

+	    m_size = _vector.size();

+	    m_vector_view = gsl_vector_view_array(&(_vector[0]), _vector.size());

+	    m_vector = &(m_vector_view.vector);

+	}

+

+	VectorView & operator=(VectorView const& _vector)

+	{

+		assert( size() == _vector.size() );

+		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());

+		return (*this);

+	}

+

+	VectorView & operator=(base_type::base_type const& _vector)

+	{

+		assert( size() == _vector.size() );

+		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());

+		return (*this);

+	}

+

+	VectorView & scale(double x)

+	{

+		return static_cast<VectorView&>( base_type::scale(x) );

+	}

+

+	VectorView & translate(double x)

+	{

+		return static_cast<VectorView&>( base_type::translate(x) );

+	}

+

+	VectorView & operator+=(base_type::base_type const& _vector)

+	{

+		return static_cast<VectorView&>( base_type::operator+=(_vector) );

+	}

+

+	VectorView & operator-=(base_type::base_type const& _vector)

+	{

+		return static_cast<VectorView&>( base_type::operator-=(_vector) );

+	}

+

+	friend

+	void swap_view(VectorView & v1, VectorView & v2);

+

+  private:

+	gsl_vector_view m_vector_view;

+

+}; // end class VectorView

+

+

+inline

+void swap_view(VectorView & v1, VectorView & v2)

+{

+	assert( v1.size() == v2.size() );

+	std::swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too

+}

+

+

+} } // end namespaces

+

+

+#endif /*_NL_VECTOR_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:encoding=utf-8:textwidth=99 :