EOL fixup

(bzr r6149)

5 files changed, 2213 insertions, 2213 deletions

diff --git a/src/2geom/numeric/fitting-model.h b/src/2geom/numeric/fitting-model.h
index 145be40e4..cc3113372 100644
--- a/src/2geom/numeric/fitting-model.h
+++ b/src/2geom/numeric/fitting-model.h
@@ -1,423 +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 :

+/*
+ * 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
index edacc663a..d589d86e3 100644
--- a/src/2geom/numeric/fitting-tool.h
+++ b/src/2geom/numeric/fitting-tool.h
@@ -1,532 +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 :

+/*
+ * 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 5b516c9e6..dc2a1d7e0 100644
--- a/src/2geom/numeric/linear_system.h
+++ b/src/2geom/numeric/linear_system.h
@@ -1,138 +1,138 @@
-/*

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

+/*
+ * 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.h b/src/2geom/numeric/matrix.h
index 64557a6f1..156b6e9a2 100644
--- a/src/2geom/numeric/matrix.h
+++ b/src/2geom/numeric/matrix.h
@@ -1,555 +1,555 @@
-/*

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

+/*
+ * 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 43a39a1ac..3e53405f4 100644
--- a/src/2geom/numeric/vector.h
+++ b/src/2geom/numeric/vector.h
@@ -1,565 +1,565 @@
-/*

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

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