summaryrefslogtreecommitdiffstats
path: root/src/2geom/numeric
diff options
context:
space:
mode:
Diffstat (limited to 'src/2geom/numeric')
-rw-r--r--src/2geom/numeric/fitting-model.h521
-rw-r--r--src/2geom/numeric/fitting-tool.h562
-rw-r--r--src/2geom/numeric/linear_system.h138
-rw-r--r--src/2geom/numeric/matrix.cpp154
-rw-r--r--src/2geom/numeric/matrix.h603
-rw-r--r--src/2geom/numeric/symmetric-matrix-fs-operation.h102
-rw-r--r--src/2geom/numeric/symmetric-matrix-fs-trace.h427
-rw-r--r--src/2geom/numeric/symmetric-matrix-fs.h733
-rw-r--r--src/2geom/numeric/vector.h594
9 files changed, 3834 insertions, 0 deletions
diff --git a/src/2geom/numeric/fitting-model.h b/src/2geom/numeric/fitting-model.h
new file mode 100644
index 0000000..0316f57
--- /dev/null
+++ b/src/2geom/numeric/fitting-model.h
@@ -0,0 +1,521 @@
+/*
+ * 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/polynomial.h>
+#include <2geom/ellipse.h>
+#include <2geom/circle.h>
+#include <2geom/utils.h>
+#include <2geom/conicsec.h>
+
+
+namespace Geom { namespace NL {
+
+/*
+ * A model is an abstraction for an expression dependent from a parameter where
+ * the coefficients of this expression are the unknowns of the fitting problem.
+ * For a ceratain number of parameter values we know the related values
+ * the expression evaluates to: from each parameter value we get a row of
+ * the matrix of the fitting problem, from each expression value we get
+ * the related constant term.
+ * Example: given the model a*x^2 + b*x + c = 0; from x = 1 we get
+ * the equation a + b + c = 0, in this example the constant term is always
+ * the same for each parameter value.
+ *
+ * A model is required to implement 3 methods:
+ *
+ * - size : returns the number of unknown coefficients that appear in
+ * the expression of the fitting problem;
+ * - feed : its input is a parameter value and the related expression value,
+ * it generates a matrix row and a new entry of the constant vector
+ * of the fitting problem;
+ * - instance : it has an input parameter represented by the raw vector
+ * solution of the fitting problem and an output parameter
+ * of type InstanceType that return a specific object that is
+ * generated using the fitting problem solution, in the example
+ * above the object could be a Poly type.
+ */
+
+/*
+ * 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
+// represented 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
+// represented 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;
+// }
+//};
+
+// 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 a conic section
+template< typename InstanceType >
+class LFMConicEquation
+ : 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 LFMConicSection
+ : public LFMConicEquation<xAx>
+{
+ public:
+ void instance(xAx & c, ConstVectorView const& coeff) const
+ {
+ c.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);
+ }
+};
+
+// this model generates Ellipse curves
+class LFMEllipse
+ : public LFMConicEquation<Ellipse>
+{
+ public:
+ void instance(Ellipse & e, ConstVectorView const& coeff) const
+ {
+ e.setCoefficients(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);
+ }
+};
+
+
+// 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 a circle curve
+template< typename InstanceType >
+class LFMCircleEquation
+ : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>
+{
+ public:
+ void feed( VectorView & coeff, double & fixed_term, Point const& p ) const
+ {
+ coeff[0] = p[X];
+ coeff[1] = p[Y];
+ coeff[2] = 1;
+ fixed_term = p[X] * p[X] + p[Y] * p[Y];
+ }
+
+ size_t size() const
+ {
+ return 3;
+ }
+};
+
+
+// this model generates Ellipse curves
+class LFMCircle
+ : public LFMCircleEquation<Circle>
+{
+ public:
+ void instance(Circle & c, ConstVectorView const& coeff) const
+ {
+ c.setCoefficients(1, coeff[0], coeff[1], coeff[2]);
+ }
+};
+
+
+// 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.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 degree>
+class LFMBezierCurveN
+ : public LinearFittingModel< double, Point, BezierCurveN<degree> >
+{
+ public:
+ LFMBezierCurveN()
+ : mob(degree+1)
+ {
+ }
+
+ void feed( VectorView & coeff, double t ) const
+ {
+ mob.feed(coeff, t);
+ }
+
+ size_t size() const
+ {
+ return mob.size();
+ }
+
+ void instance(BezierCurveN<degree> & bc, ConstMatrixView const& raw_data) const
+ {
+ Bezier bx(degree);
+ Bezier by(degree);
+ mob.instance(bx, raw_data.column_const_view(X));
+ mob.instance(by, raw_data.column_const_view(Y));
+ bc = BezierCurveN<degree>(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:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/fitting-tool.h b/src/2geom/numeric/fitting-tool.h
new file mode 100644
index 0000000..f2e856a
--- /dev/null
+++ b/src/2geom/numeric/fitting-tool.h
@@ -0,0 +1,562 @@
+/*
+ * Fitting Tools
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _NL_FITTING_TOOL_H_
+#define _NL_FITTING_TOOL_H_
+
+
+#include <2geom/numeric/vector.h>
+#include <2geom/numeric/matrix.h>
+
+#include <2geom/point.h>
+
+#include <vector>
+
+
+/*
+ * The least_square_fitter class represents a tool for solving a fitting
+ * problem with respect to a given model that represents an expression
+ * dependent from a parameter where the coefficients of this expression
+ * are the unknowns of the fitting problem.
+ * The minimizing solution is found by computing the pseudo-inverse
+ * of the problem matrix
+ */
+
+
+namespace Geom { namespace NL {
+
+namespace detail {
+
+template< typename ModelT>
+class lsf_base
+{
+ public:
+ typedef ModelT model_type;
+ typedef typename model_type::parameter_type parameter_type;
+ typedef typename model_type::value_type value_type;
+
+ lsf_base( model_type const& _model, size_t forecasted_samples )
+ : m_model(_model),
+ m_total_samples(0),
+ m_matrix(forecasted_samples, m_model.size()),
+ m_psdinv_matrix(NULL)
+ {
+ }
+
+ // compute pseudo inverse
+ void update()
+ {
+ if (total_samples() == 0) return;
+ if (m_psdinv_matrix != NULL)
+ {
+ delete m_psdinv_matrix;
+ }
+ MatrixView mv(m_matrix, 0, 0, total_samples(), m_matrix.columns());
+ m_psdinv_matrix = new Matrix( pseudo_inverse(mv) );
+ assert(m_psdinv_matrix != NULL);
+ }
+
+ size_t total_samples() const
+ {
+ return m_total_samples;
+ }
+
+ bool is_full() const
+ {
+ return (total_samples() == m_matrix.rows());
+ }
+
+ void clear()
+ {
+ m_total_samples = 0;
+ }
+
+ virtual
+ ~lsf_base()
+ {
+ if (m_psdinv_matrix != NULL)
+ {
+ delete m_psdinv_matrix;
+ }
+ }
+
+ protected:
+ const model_type & m_model;
+ size_t m_total_samples;
+ Matrix m_matrix;
+ Matrix* m_psdinv_matrix;
+
+}; // end class lsf_base
+
+
+
+
+template< typename ModelT, typename ValueType = typename ModelT::value_type>
+class lsf_solution
+{
+};
+
+// a fitting process on samples with value of type double
+// produces a solution of type Vector
+template< typename ModelT>
+class lsf_solution<ModelT, double>
+ : public lsf_base<ModelT>
+{
+public:
+ typedef ModelT model_type;
+ typedef typename model_type::parameter_type parameter_type;
+ typedef typename model_type::value_type value_type;
+ typedef Vector solution_type;
+ typedef lsf_base<model_type> base_type;
+
+ using base_type::m_model;
+ using base_type::m_psdinv_matrix;
+ using base_type::total_samples;
+
+public:
+ lsf_solution<ModelT, double>( model_type const& _model,
+ size_t forecasted_samples )
+ : base_type(_model, forecasted_samples),
+ m_solution(_model.size())
+ {
+ }
+
+ template< typename VectorT >
+ solution_type& result(VectorT const& sample_values)
+ {
+ assert(sample_values.size() == total_samples());
+ ConstVectorView sv(sample_values);
+ m_solution = (*m_psdinv_matrix) * sv;
+ return m_solution;
+ }
+
+ // a comparison between old sample values and the new ones is performed
+ // in order to minimize computation
+ // prerequisite:
+ // old_sample_values.size() == new_sample_values.size()
+ // no update() call can be performed between two result invocations
+ template< typename VectorT >
+ solution_type& result( VectorT const& old_sample_values,
+ VectorT const& new_sample_values )
+ {
+ assert(old_sample_values.size() == total_samples());
+ assert(new_sample_values.size() == total_samples());
+ Vector diff(total_samples());
+ for (size_t i = 0; i < diff.size(); ++i)
+ {
+ diff[i] = new_sample_values[i] - old_sample_values[i];
+ }
+ Vector column(m_model.size());
+ Vector delta(m_model.size(), 0.0);
+ for (size_t i = 0; i < diff.size(); ++i)
+ {
+ if (diff[i] != 0)
+ {
+ column = m_psdinv_matrix->column_view(i);
+ column.scale(diff[i]);
+ delta += column;
+ }
+ }
+ m_solution += delta;
+ return m_solution;
+ }
+
+ solution_type& result()
+ {
+ return m_solution;
+ }
+
+private:
+ solution_type m_solution;
+
+}; // end class lsf_solution<ModelT, double>
+
+
+// a fitting process on samples with value of type Point
+// produces a solution of type Matrix (with 2 columns)
+template< typename ModelT>
+class lsf_solution<ModelT, Point>
+ : public lsf_base<ModelT>
+{
+public:
+ typedef ModelT model_type;
+ typedef typename model_type::parameter_type parameter_type;
+ typedef typename model_type::value_type value_type;
+ typedef Matrix solution_type;
+ typedef lsf_base<model_type> base_type;
+
+ using base_type::m_model;
+ using base_type::m_psdinv_matrix;
+ using base_type::total_samples;
+
+public:
+ lsf_solution<ModelT, Point>( model_type const& _model,
+ size_t forecasted_samples )
+ : base_type(_model, forecasted_samples),
+ m_solution(_model.size(), 2)
+ {
+ }
+
+ solution_type& result(std::vector<Point> const& sample_values)
+ {
+ assert(sample_values.size() == total_samples());
+ Matrix svm(total_samples(), 2);
+ for (size_t i = 0; i < total_samples(); ++i)
+ {
+ svm(i, X) = sample_values[i][X];
+ svm(i, Y) = sample_values[i][Y];
+ }
+ m_solution = (*m_psdinv_matrix) * svm;
+ return m_solution;
+ }
+
+ // a comparison between old sample values and the new ones is performed
+ // in order to minimize computation
+ // prerequisite:
+ // old_sample_values.size() == new_sample_values.size()
+ // no update() call can to be performed between two result invocations
+ solution_type& result( std::vector<Point> const& old_sample_values,
+ std::vector<Point> const& new_sample_values )
+ {
+ assert(old_sample_values.size() == total_samples());
+ assert(new_sample_values.size() == total_samples());
+ Matrix diff(total_samples(), 2);
+ for (size_t i = 0; i < total_samples(); ++i)
+ {
+ diff(i, X) = new_sample_values[i][X] - old_sample_values[i][X];
+ diff(i, Y) = new_sample_values[i][Y] - old_sample_values[i][Y];
+ }
+ Vector column(m_model.size());
+ Matrix delta(m_model.size(), 2, 0.0);
+ VectorView deltax = delta.column_view(X);
+ VectorView deltay = delta.column_view(Y);
+ for (size_t i = 0; i < total_samples(); ++i)
+ {
+ if (diff(i, X) != 0)
+ {
+ column = m_psdinv_matrix->column_view(i);
+ column.scale(diff(i, X));
+ deltax += column;
+ }
+ if (diff(i, Y) != 0)
+ {
+ column = m_psdinv_matrix->column_view(i);
+ column.scale(diff(i, Y));
+ deltay += column;
+ }
+ }
+ m_solution += delta;
+ return m_solution;
+ }
+
+ solution_type& result()
+ {
+ return m_solution;
+ }
+
+private:
+ solution_type m_solution;
+
+}; // end class lsf_solution<ModelT, Point>
+
+
+
+
+template< typename ModelT,
+ bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
+class lsf_with_fixed_terms
+{
+};
+
+
+// fitting tool for completely unknown models
+template< typename ModelT>
+class lsf_with_fixed_terms<ModelT, false>
+ : public lsf_solution<ModelT>
+{
+ public:
+ typedef ModelT model_type;
+ typedef typename model_type::parameter_type parameter_type;
+ typedef typename model_type::value_type value_type;
+ typedef lsf_solution<model_type> base_type;
+ typedef typename base_type::solution_type solution_type;
+
+ using base_type::total_samples;
+ using base_type::is_full;
+ using base_type::m_matrix;
+ using base_type::m_total_samples;
+ using base_type::m_model;
+
+ public:
+ lsf_with_fixed_terms<ModelT, false>( model_type const& _model,
+ size_t forecasted_samples )
+ : base_type(_model, forecasted_samples)
+ {
+ }
+
+ void append(parameter_type const& sample_parameter)
+ {
+ assert(!is_full());
+ VectorView row = m_matrix.row_view(total_samples());
+ m_model.feed(row, sample_parameter);
+ ++m_total_samples;
+ }
+
+ void append_copy(size_t sample_index)
+ {
+ assert(!is_full());
+ assert(sample_index < total_samples());
+ VectorView dest_row = m_matrix.row_view(total_samples());
+ VectorView source_row = m_matrix.row_view(sample_index);
+ dest_row = source_row;
+ ++m_total_samples;
+ }
+
+}; // end class lsf_with_fixed_terms<ModelT, false>
+
+
+// fitting tool for partially known models
+template< typename ModelT>
+class lsf_with_fixed_terms<ModelT, true>
+ : public lsf_solution<ModelT>
+{
+ public:
+ typedef ModelT model_type;
+ typedef typename model_type::parameter_type parameter_type;
+ typedef typename model_type::value_type value_type;
+ typedef lsf_solution<model_type> base_type;
+ typedef typename base_type::solution_type solution_type;
+
+ using base_type::total_samples;
+ using base_type::is_full;
+ using base_type::m_matrix;
+ using base_type::m_total_samples;
+ using base_type::m_model;
+
+ public:
+ lsf_with_fixed_terms<ModelT, true>( model_type const& _model,
+ size_t forecasted_samples )
+ : base_type(_model, forecasted_samples),
+ m_vector(forecasted_samples),
+ m_vector_view(NULL)
+ {
+ }
+ void append(parameter_type const& sample_parameter)
+ {
+ assert(!is_full());
+ VectorView row = m_matrix.row_view(total_samples());
+ m_model.feed(row, m_vector[total_samples()], sample_parameter);
+ ++m_total_samples;
+ }
+
+ void append_copy(size_t sample_index)
+ {
+ assert(!is_full());
+ assert(sample_index < total_samples());
+ VectorView dest_row = m_matrix.row_view(total_samples());
+ VectorView source_row = m_matrix.row_view(sample_index);
+ dest_row = source_row;
+ m_vector[total_samples()] = m_vector[sample_index];
+ ++m_total_samples;
+ }
+
+ void update()
+ {
+ base_type::update();
+ if (total_samples() == 0) return;
+ if (m_vector_view != NULL)
+ {
+ delete m_vector_view;
+ }
+ m_vector_view = new VectorView(m_vector, base_type::total_samples());
+ assert(m_vector_view != NULL);
+ }
+
+ 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; // VSC legacy support
+ solution_type& result( std::vector<Point> const& old_sample_values,
+ std::vector<Point> const& new_sample_values )
+ {
+ return base_type::result(old_sample_values, new_sample_values);
+ }
+
+ solution_type& result()
+ {
+ return base_type::result();
+ }
+
+ public:
+ least_squeares_fitter<ModelT, double, true>( model_type const& _model,
+ size_t forecasted_samples )
+ : base_type(_model, forecasted_samples)
+ {
+ }
+
+ template< typename VectorT >
+ solution_type& result(VectorT const& sample_values)
+ {
+ assert(sample_values.size() == m_vector_view->size());
+ Vector sv(sample_values.size());
+ for (size_t i = 0; i < sv.size(); ++i)
+ sv[i] = sample_values[i] - (*m_vector_view)[i];
+ return base_type::result(sv);
+ }
+
+}; // end class least_squeares_fitter<ModelT, double, true>
+
+
+template< typename ModelT>
+class least_squeares_fitter<ModelT, Point, true>
+ : public detail::lsf_with_fixed_terms<ModelT>
+{
+ public:
+ typedef ModelT model_type;
+ typedef detail::lsf_with_fixed_terms<model_type> base_type;
+ typedef typename base_type::parameter_type parameter_type;
+ typedef typename base_type::value_type value_type;
+ typedef typename base_type::solution_type solution_type;
+
+ using base_type::m_vector_view;
+ //using base_type::result; // VCS legacy support
+ solution_type& result( std::vector<Point> const& old_sample_values,
+ std::vector<Point> const& new_sample_values )
+ {
+ return base_type::result(old_sample_values, new_sample_values);
+ }
+
+ solution_type& result()
+ {
+ return base_type::result();
+ }
+
+
+ public:
+ least_squeares_fitter<ModelT, Point, true>( model_type const& _model,
+ size_t forecasted_samples )
+ : base_type(_model, forecasted_samples)
+ {
+ }
+
+ solution_type& result(std::vector<Point> const& sample_values)
+ {
+ assert(sample_values.size() == m_vector_view->size());
+ NL::Matrix sv(sample_values.size(), 2);
+ for (size_t i = 0; i < sample_values.size(); ++i)
+ {
+ sv(i, X) = sample_values[i][X] - (*m_vector_view)[i];
+ sv(i, Y) = sample_values[i][Y] - (*m_vector_view)[i];
+ }
+ return base_type::result(sv);
+ }
+
+}; // end class least_squeares_fitter<ModelT, Point, true>
+
+
+} // end namespace NL
+} // end namespace Geom
+
+
+
+#endif // _NL_FITTING_TOOL_H_
+
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/linear_system.h b/src/2geom/numeric/linear_system.h
new file mode 100644
index 0000000..f793e20
--- /dev/null
+++ b/src/2geom/numeric/linear_system.h
@@ -0,0 +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:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/matrix.cpp b/src/2geom/numeric/matrix.cpp
new file mode 100644
index 0000000..98ff3b6
--- /dev/null
+++ b/src/2geom/numeric/matrix.cpp
@@ -0,0 +1,154 @@
+/*
+ * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines;
+ * "views" mimic the semantic of C++ references: any operation performed
+ * on a "view" is actually performed on the "viewed object"
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#include <2geom/numeric/matrix.h>
+#include <2geom/numeric/vector.h>
+
+
+namespace Geom { namespace NL {
+
+Vector operator*( detail::BaseMatrixImpl const& A,
+ detail::BaseVectorImpl const& v )
+{
+ assert(A.columns() == v.size());
+
+ Vector result(A.rows(), 0.0);
+ for (size_t i = 0; i < A.rows(); ++i)
+ for (size_t j = 0; j < A.columns(); ++j)
+ result[i] += A(i,j) * v[j];
+
+ return result;
+}
+
+Matrix operator*( detail::BaseMatrixImpl const& A,
+ detail::BaseMatrixImpl const& B )
+{
+ assert(A.columns() == B.rows());
+
+ Matrix C(A.rows(), B.columns(), 0.0);
+ for (size_t i = 0; i < C.rows(); ++i)
+ for (size_t j = 0; j < C.columns(); ++j)
+ for (size_t k = 0; k < A.columns(); ++k)
+ C(i,j) += A(i,k) * B(k, j);
+
+ return C;
+}
+
+Matrix pseudo_inverse(detail::BaseMatrixImpl const& A)
+{
+
+ Matrix U(A);
+ Matrix V(A.columns(), A.columns());
+ Vector s(A.columns());
+ gsl_vector* work = gsl_vector_alloc(A.columns());
+
+ gsl_linalg_SV_decomp( U.get_gsl_matrix(),
+ V.get_gsl_matrix(),
+ s.get_gsl_vector(),
+ work );
+
+ Matrix P(A.columns(), A.rows(), 0.0);
+
+ int sz = s.size();
+ while ( sz-- > 0 && s[sz] == 0 ) {}
+ ++sz;
+ if (sz == 0) return P;
+ VectorView sv(s, sz);
+
+ for (size_t i = 0; i < sv.size(); ++i)
+ {
+ VectorView v = V.column_view(i);
+ v.scale(1/sv[i]);
+ for (size_t h = 0; h < P.rows(); ++h)
+ for (size_t k = 0; k < P.columns(); ++k)
+ P(h,k) += V(h,i) * U(k,i);
+ }
+
+ return P;
+}
+
+
+double trace (detail::BaseMatrixImpl const& A)
+{
+ if (A.rows() != A.columns())
+ {
+ THROW_RANGEERROR ("NL::Matrix: computing trace: "
+ "rows() != columns()");
+ }
+ double t = 0;
+ for (size_t i = 0; i < A.rows(); ++i)
+ {
+ t += A(i,i);
+ }
+ return t;
+}
+
+
+double det (detail::BaseMatrixImpl const& A)
+{
+ if (A.rows() != A.columns())
+ {
+ THROW_RANGEERROR ("NL::Matrix: computing determinant: "
+ "rows() != columns()");
+ }
+
+ Matrix LU(A);
+ int s;
+ gsl_permutation * p = gsl_permutation_alloc(LU.rows());
+ gsl_linalg_LU_decomp (LU.get_gsl_matrix(), p, &s);
+
+ double t = 1;
+ for (size_t i = 0; i < LU.rows(); ++i)
+ {
+ t *= LU(i,i);
+ }
+
+ gsl_permutation_free(p);
+ return t;
+}
+
+
+} } // end namespaces
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/matrix.h b/src/2geom/numeric/matrix.h
new file mode 100644
index 0000000..5cebf10
--- /dev/null
+++ b/src/2geom/numeric/matrix.h
@@ -0,0 +1,603 @@
+/*
+ * 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/exception.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
+ {
+ 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;
+ }
+
+ bool is_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;
+ }
+
+ 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(); // VSC legacy support
+ const double & operator() (size_t i, size_t j) const
+ {
+ return base_type::operator ()(i, j);
+ }
+
+ 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<<;
+
+
+template <size_t N>
+class ConstBaseSymmetricMatrix;
+
+
+class Matrix: public detail::MatrixImpl
+{
+ public:
+ typedef detail::MatrixImpl base_type;
+
+ public:
+ // the matrix is not initialized
+ 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());
+ }
+
+ template <size_t N>
+ explicit
+ Matrix(ConstBaseSymmetricMatrix<N> const& _smatrix)
+ {
+ m_rows = N;
+ m_columns = N;
+ m_matrix = gsl_matrix_alloc(N, N);
+ for (size_t i = 0; i < N; ++i)
+ for (size_t j = 0; j < N ; ++j)
+ (*gsl_matrix_ptr(m_matrix, i, j)) = _smatrix(i,j);
+ }
+
+ 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;
+ }
+
+ template <size_t N>
+ Matrix & operator=(ConstBaseSymmetricMatrix<N> const& _smatrix)
+ {
+ assert (rows() == N && columns() == N);
+ for (size_t i = 0; i < N; ++i)
+ for (size_t j = 0; j < N ; ++j)
+ (*this)(i,j) = _smatrix(i,j);
+ 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());
+ using std::swap;
+ swap(m1.m_matrix, m2.m_matrix);
+}
+
+inline void swap_any(Matrix &m1, Matrix &m2)
+{
+ using std::swap;
+ swap(m1.m_matrix, m2.m_matrix);
+ swap(m1.m_rows, m2.m_rows);
+ 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());
+ using std::swap;
+ 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);
+
+double trace (detail::BaseMatrixImpl const& A);
+
+double det (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:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/symmetric-matrix-fs-operation.h b/src/2geom/numeric/symmetric-matrix-fs-operation.h
new file mode 100644
index 0000000..c5aaa72
--- /dev/null
+++ b/src/2geom/numeric/symmetric-matrix-fs-operation.h
@@ -0,0 +1,102 @@
+/*
+ * SymmetricMatrix basic operation
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2009 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_SYMMETRIC_MATRIX_FS_OPERATION_H_
+#define _NL_SYMMETRIC_MATRIX_FS_OPERATION_H_
+
+
+#include <2geom/numeric/symmetric-matrix-fs.h>
+#include <2geom/numeric/symmetric-matrix-fs-trace.h>
+
+
+
+
+namespace Geom { namespace NL {
+
+template <size_t N>
+SymmetricMatrix<N> adj(const ConstBaseSymmetricMatrix<N> & S);
+
+template <>
+inline
+SymmetricMatrix<2> adj(const ConstBaseSymmetricMatrix<2> & S)
+{
+ SymmetricMatrix<2> result;
+ result.get<0,0>() = S.get<1,1>();
+ result.get<1,0>() = -S.get<1,0>();
+ result.get<1,1>() = S.get<0,0>();
+ return result;
+}
+
+template <>
+inline
+SymmetricMatrix<3> adj(const ConstBaseSymmetricMatrix<3> & S)
+{
+ SymmetricMatrix<3> result;
+
+ result.get<0,0>() = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>();
+ result.get<1,0>() = S.get<0,2>() * S.get<2,1>() - S.get<0,1>() * S.get<2,2>();
+ result.get<1,1>() = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>();
+ result.get<2,0>() = S.get<0,1>() * S.get<1,2>() - S.get<0,2>() * S.get<1,1>();
+ result.get<2,1>() = S.get<0,2>() * S.get<1,0>() - S.get<0,0>() * S.get<1,2>();
+ result.get<2,2>() = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
+ return result;
+}
+
+template <size_t N>
+inline
+SymmetricMatrix<N> inverse(const ConstBaseSymmetricMatrix<N> & S)
+{
+ SymmetricMatrix<N> result = adj(S);
+ double d = det(S);
+ assert (d != 0);
+ result.scale (1/d);
+ return result;
+}
+
+} /* end namespace NL*/ } /* end namespace Geom*/
+
+
+#endif // _NL_SYMMETRIC_MATRIX_FS_OPERATION_H_
+
+
+
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/symmetric-matrix-fs-trace.h b/src/2geom/numeric/symmetric-matrix-fs-trace.h
new file mode 100644
index 0000000..eff3dd2
--- /dev/null
+++ b/src/2geom/numeric/symmetric-matrix-fs-trace.h
@@ -0,0 +1,427 @@
+/*
+ * SymmetricMatrix trace
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2009 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_TRACE_H_
+#define _NL_TRACE_H_
+
+
+#include <2geom/numeric/matrix.h>
+#include <2geom/numeric/symmetric-matrix-fs.h>
+
+
+
+
+
+namespace Geom { namespace NL {
+
+
+namespace detail
+{
+
+/*
+ * helper routines
+ */
+
+inline
+int sgn_prod (int x, int y)
+{
+ if (x == 0 || y == 0) return 0;
+ if (x == y) return 1;
+ return -1;
+}
+
+inline
+bool abs_less (double x, double y)
+{
+ return (std::fabs(x) < std::fabs(y));
+}
+
+
+/*
+ * trace K-th of symmetric matrix S of order N
+ */
+template <size_t K, size_t N>
+struct trace
+{
+ static double evaluate(const ConstBaseSymmetricMatrix<N> &S);
+};
+
+template <size_t N>
+struct trace<1,N>
+{
+ static
+ double evaluate (const ConstBaseSymmetricMatrix<N> & S)
+ {
+ double t = 0;
+ for (size_t i = 0; i < N; ++i)
+ {
+ t += S(i,i);
+ }
+ return t;
+ }
+};
+
+template <size_t N>
+struct trace<N,N>
+{
+ static
+ double evaluate (const ConstBaseSymmetricMatrix<N> & S)
+ {
+ Matrix M(S);
+ return det(M);
+ }
+};
+
+/*
+ * trace for symmetric matrix of order 2
+ */
+template <>
+struct trace<1,2>
+{
+ static
+ double evaluate (const ConstBaseSymmetricMatrix<2> & S)
+ {
+ return (S.get<0,0>() + S.get<1,1>());
+ }
+};
+
+template <>
+struct trace<2,2>
+{
+ static
+ double evaluate (const ConstBaseSymmetricMatrix<2> & S)
+ {
+ return (S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>());
+ }
+};
+
+
+/*
+ * trace for symmetric matrix of order 3
+ */
+template <>
+struct trace<1,3>
+{
+ static
+ double evaluate (const ConstBaseSymmetricMatrix<3> & S)
+ {
+ return (S.get<0,0>() + S.get<1,1>() + S.get<2,2>());
+ }
+};
+
+template <>
+struct trace<2,3>
+{
+ static
+ double evaluate (const ConstBaseSymmetricMatrix<3> & S)
+ {
+ double a00 = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>();
+ double a11 = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>();
+ double a22 = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
+ return (a00 + a11 + a22);
+ }
+};
+
+template <>
+struct trace<3,3>
+{
+ static
+ double evaluate (const ConstBaseSymmetricMatrix<3> & S)
+ {
+ double d = S.get<0,0>() * S.get<1,1>() * S.get<2,2>();
+ d += (2 * S.get<1,0>() * S.get<2,0>() * S.get<2,1>());
+ d -= (S.get<0,0>() * S.get<2,1>() * S.get<2,1>());
+ d -= (S.get<1,1>() * S.get<2,0>() * S.get<2,0>());
+ d -= (S.get<2,2>() * S.get<1,0>() * S.get<1,0>());
+ return d;
+ }
+};
+
+
+/*
+ * sign of trace K-th
+ */
+template <size_t K, size_t N>
+struct trace_sgn
+{
+ static
+ int evaluate (const ConstBaseSymmetricMatrix<N> & S)
+ {
+ double d = trace<K, N>::evaluate(S);
+ return sgn(d);
+ }
+};
+
+
+/*
+ * sign of trace for symmetric matrix of order 2
+ */
+template <>
+struct trace_sgn<2,2>
+{
+ static
+ int evaluate (const ConstBaseSymmetricMatrix<2> & S)
+ {
+ double m00 = S.get<0,0>();
+ double m10 = S.get<1,0>();
+ double m11 = S.get<1,1>();
+
+ int sm00 = sgn (m00);
+ int sm10 = sgn (m10);
+ int sm11 = sgn (m11);
+
+ if (sm10 == 0)
+ {
+ return sgn_prod (sm00, sm11);
+ }
+ else
+ {
+ int sm00m11 = sgn_prod (sm00, sm11);
+ if (sm00m11 == 1)
+ {
+ int e00, e10, e11;
+ double f00 = std::frexp (m00, &e00);
+ double f10 = std::frexp (m10, &e10);
+ double f11 = std::frexp (m11, &e11);
+
+ int e0011 = e00 + e11;
+ int e1010 = e10 << 1;
+ int ed = e0011 - e1010;
+
+ if (ed > 1)
+ {
+ return 1;
+ }
+ else if (ed < -1)
+ {
+ return -1;
+ }
+ else
+ {
+ double d = std::ldexp (f00 * f11, ed) - f10 * f10;
+ //std::cout << "trace_sgn<2,2>: det = " << d << std::endl;
+ double eps = std::ldexp (1, -50);
+ if (std::fabs(d) < eps) return 0;
+ return sgn (d);
+ }
+ }
+ return -1;
+ }
+ }
+};
+
+
+/*
+ * sign of trace for symmetric matrix of order 3
+ */
+template <>
+struct trace_sgn<2,3>
+{
+ static
+ int evaluate (const ConstBaseSymmetricMatrix<3> & S)
+ {
+ double eps = std::ldexp (1, -50);
+ double t[6];
+
+ t[0] = S.get<1,1>() * S.get<2,2>();
+ t[1] = - S.get<1,2>() * S.get<2,1>();
+ t[2] = S.get<0,0>() * S.get<2,2>();
+ t[3] = - S.get<0,2>() * S.get<2,0>();
+ t[4] = S.get<0,0>() * S.get<1,1>();
+ t[5] = - S.get<0,1>() * S.get<1,0>();
+
+
+ double* maxp = std::max_element (t, t+6, abs_less);
+ int em;
+ std::frexp(*maxp, &em);
+ double d = 0;
+ for (size_t i = 0; i < 6; ++i)
+ {
+ d += t[i];
+ }
+ double r = std::fabs (std::ldexp (d, -em)); // relative error
+ //std::cout << "trace_sgn<2,3>: d = " << d << std::endl;
+ //std::cout << "trace_sgn<2,3>: r = " << r << std::endl;
+ if (r < eps) return 0;
+ if (d > 0) return 1;
+ return -1;
+ }
+};
+
+template <>
+struct trace_sgn<3,3>
+{
+ static
+ int evaluate (const ConstBaseSymmetricMatrix<3> & S)
+ {
+
+ double eps = std::ldexp (1, -48);
+ double t[5];
+
+ t[0] = S.get<0,0>() * S.get<1,1>() * S.get<2,2>();
+ t[1] = 2 * S.get<1,0>() * S.get<2,0>() * S.get<2,1>();
+ t[2] = -(S.get<0,0>() * S.get<2,1>() * S.get<2,1>());
+ t[3] = -(S.get<1,1>() * S.get<2,0>() * S.get<2,0>());
+ t[4] = -(S.get<2,2>() * S.get<1,0>() * S.get<1,0>());
+
+ double* maxp = std::max_element (t, t+5, abs_less);
+ int em;
+ std::frexp(*maxp, &em);
+ double d = 0;
+ for (size_t i = 0; i < 5; ++i)
+ {
+ d += t[i];
+ }
+ //std::cout << "trace_sgn<3,3>: d = " << d << std::endl;
+ double r = std::fabs (std::ldexp (d, -em)); // relative error
+ //std::cout << "trace_sgn<3,3>: r = " << r << std::endl;
+
+ if (r < eps) return 0;
+ if (d > 0) return 1;
+ return -1;
+ }
+}; // end struct trace_sgn<3,3>
+
+} // end namespace detail
+
+
+template <size_t K, size_t N>
+inline
+double trace (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+ return detail::trace<K, N>::evaluate(_matrix);
+}
+
+template <size_t N>
+inline
+double trace (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+ return detail::trace<1, N>::evaluate(_matrix);
+}
+
+template <size_t N>
+inline
+double det (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+ return detail::trace<N, N>::evaluate(_matrix);
+}
+
+
+template <size_t K, size_t N>
+inline
+int trace_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+ return detail::trace_sgn<K, N>::evaluate(_matrix);
+}
+
+template <size_t N>
+inline
+int trace_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+ return detail::trace_sgn<1, N>::evaluate(_matrix);
+}
+
+template <size_t N>
+inline
+int det_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+ return detail::trace_sgn<N, N>::evaluate(_matrix);
+}
+
+/*
+template <size_t N>
+inline
+size_t rank (const ConstBaseSymmetricMatrix<N> & S)
+{
+ THROW_NOTIMPLEMENTED();
+ return 0;
+}
+
+template <>
+inline
+size_t rank<2> (const ConstBaseSymmetricMatrix<2> & S)
+{
+ if (S.is_zero()) return 0;
+ double d = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
+ if (d != 0) return 2;
+ return 1;
+}
+
+template <>
+inline
+size_t rank<3> (const ConstBaseSymmetricMatrix<3> & S)
+{
+ if (S.is_zero()) return 0;
+
+ double a20 = S.get<0,1>() * S.get<1,2>() - S.get<0,2>() * S.get<1,1>();
+ double a21 = S.get<0,2>() * S.get<1,0>() - S.get<0,0>() * S.get<1,2>();
+ double a22 = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
+ double d = a20 * S.get<2,0>() + a21 * S.get<2,1>() + a22 * S.get<2,2>();
+
+ if (d != 0) return 3;
+
+ if (a20 != 0 || a21 != 0 || a22 != 0) return 2;
+
+ double a00 = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>();
+ if (a00 != 0) return 2;
+
+ double a10 = S.get<0,2>() * S.get<2,1>() - S.get<0,1>() * S.get<2,2>();
+ if (a10 != 0) return 2;
+
+ double a11 = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>();
+ if (a11 != 0) return 2;
+
+ return 1;
+}
+*/
+
+} /* end namespace NL*/ } /* end namespace Geom*/
+
+
+
+
+#endif // _NL_TRACE_H_
+
+
+
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/symmetric-matrix-fs.h b/src/2geom/numeric/symmetric-matrix-fs.h
new file mode 100644
index 0000000..2fadd69
--- /dev/null
+++ b/src/2geom/numeric/symmetric-matrix-fs.h
@@ -0,0 +1,733 @@
+/*
+ * SymmetricMatrix, ConstSymmetricMatrixView, SymmetricMatrixView template
+ * classes implement fixed size symmetric matrix; "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 2009 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_SYMMETRIC_MATRIX_FS_H_
+#define _NL_SYMMETRIC_MATRIX_FS_H_
+
+
+#include <2geom/numeric/vector.h>
+#include <2geom/numeric/matrix.h>
+#include <2geom/utils.h>
+#include <2geom/exception.h>
+
+#include <boost/static_assert.hpp>
+
+#include <cassert>
+#include <utility> // for std::pair
+#include <algorithm> // for std::swap, std::copy
+#include <sstream>
+#include <string>
+
+
+
+namespace Geom { namespace NL {
+
+
+namespace detail
+{
+
+template <size_t I, size_t J, bool B = (I < J)>
+struct index
+{
+ static const size_t K = index<J, I>::K;
+};
+
+template <size_t I, size_t J>
+struct index<I, J, false>
+{
+ static const size_t K = (((I+1) * I) >> 1) + J;
+};
+
+} // end namespace detail
+
+
+
+
+template <size_t N>
+class ConstBaseSymmetricMatrix;
+
+template <size_t N>
+class BaseSymmetricMatrix;
+
+template <size_t N>
+class SymmetricMatrix;
+
+template <size_t N>
+class ConstSymmetricMatrixView;
+
+template <size_t N>
+class SymmetricMatrixView;
+
+
+
+// declaration needed for friend clause
+template <size_t N>
+bool operator== (ConstBaseSymmetricMatrix<N> const& _smatrix1,
+ ConstBaseSymmetricMatrix<N> const& _smatrix2);
+
+
+
+
+template <size_t N>
+class ConstBaseSymmetricMatrix
+{
+ public:
+ const static size_t DIM = N;
+ const static size_t DATA_SIZE = ((DIM+1) * DIM) / 2;
+
+ public:
+
+ ConstBaseSymmetricMatrix (VectorView const& _data)
+ : m_data(_data)
+ {
+ }
+
+ double operator() (size_t i, size_t j) const
+ {
+ return m_data[get_index(i,j)];
+ }
+
+ template <size_t I, size_t J>
+ double get() const
+ {
+ BOOST_STATIC_ASSERT ((I < N && J < N));
+ return m_data[detail::index<I, J>::K];
+ }
+
+
+ size_t rows() const
+ {
+ return DIM;
+ }
+
+ size_t columns() const
+ {
+ return DIM;
+ }
+
+ bool is_zero() const
+ {
+ return m_data.is_zero();
+ }
+
+ bool is_positive() const
+ {
+ return m_data.is_positive();
+ }
+
+ bool is_negative() const
+ {
+ return m_data.is_negative();
+ }
+
+ bool is_non_negative() const
+ {
+ return m_data.is_non_negative();
+ }
+
+ double min() const
+ {
+ return m_data.min();
+ }
+
+ double max() const
+ {
+ return m_data.max();
+ }
+
+ std::pair<size_t, size_t>
+ min_index() const
+ {
+ std::pair<size_t, size_t> indices(0,0);
+ double min_value = m_data[0];
+ for (size_t i = 1; i < DIM; ++i)
+ {
+ for (size_t j = 0; j <= i; ++j)
+ {
+ if (min_value > (*this)(i,j))
+ {
+ min_value = (*this)(i,j);
+ indices.first = i;
+ indices.second = j;
+ }
+ }
+ }
+ return indices;
+ }
+
+ std::pair<size_t, size_t>
+ max_index() const
+ {
+ std::pair<size_t, size_t> indices(0,0);
+ double max_value = m_data[0];
+ for (size_t i = 1; i < DIM; ++i)
+ {
+ for (size_t j = 0; j <= i; ++j)
+ {
+ if (max_value < (*this)(i,j))
+ {
+ max_value = (*this)(i,j);
+ indices.first = i;
+ indices.second = j;
+ }
+ }
+ }
+ return indices;
+ }
+
+ size_t min_on_row_index (size_t i) const
+ {
+ size_t idx = 0;
+ double min_value = (*this)(i,0);
+ for (size_t j = 1; j < DIM; ++j)
+ {
+ if (min_value > (*this)(i,j))
+ {
+ min_value = (*this)(i,j);
+ idx = j;
+ }
+ }
+ return idx;
+ }
+
+ size_t max_on_row_index (size_t i) const
+ {
+ size_t idx = 0;
+ double max_value = (*this)(i,0);
+ for (size_t j = 1; j < DIM; ++j)
+ {
+ if (max_value < (*this)(i,j))
+ {
+ max_value = (*this)(i,j);
+ idx = j;
+ }
+ }
+ return idx;
+ }
+
+ size_t min_on_column_index (size_t j) const
+ {
+ return min_on_row_index(j);
+ }
+
+ size_t max_on_column_index (size_t j) const
+ {
+ return max_on_row_index(j);
+ }
+
+ size_t min_on_diag_index () const
+ {
+ size_t idx = 0;
+ double min_value = (*this)(0,0);
+ for (size_t i = 1; i < DIM; ++i)
+ {
+ if (min_value > (*this)(i,i))
+ {
+ min_value = (*this)(i,i);
+ idx = i;
+ }
+ }
+ return idx;
+ }
+
+ size_t max_on_diag_index () const
+ {
+ size_t idx = 0;
+ double max_value = (*this)(0,0);
+ for (size_t i = 1; i < DIM; ++i)
+ {
+ if (max_value < (*this)(i,i))
+ {
+ max_value = (*this)(i,i);
+ idx = i;
+ }
+ }
+ return idx;
+ }
+
+ std::string str() const;
+
+ ConstSymmetricMatrixView<N-1> main_minor_const_view() const;
+
+ SymmetricMatrix<N> operator- () const;
+
+ Vector operator* (ConstVectorView _vector) const
+ {
+ assert (_vector.size() == DIM);
+ Vector result(DIM, 0.0);
+
+ for (size_t i = 0; i < DIM; ++i)
+ {
+ for (size_t j = 0; j < DIM; ++j)
+ {
+ result[i] += (*this)(i,j) * _vector[j];
+ }
+ }
+ return result;
+ }
+
+ protected:
+ static size_t get_index (size_t i, size_t j)
+ {
+ if (i < j) return get_index (j, i);
+ size_t k = (i+1) * i;
+ k >>= 1;
+ k += j;
+ return k;
+ }
+
+ protected:
+ ConstVectorView get_data() const
+ {
+ return m_data;
+ }
+
+ friend
+ bool operator==<N> (ConstBaseSymmetricMatrix const& _smatrix1,
+ ConstBaseSymmetricMatrix const& _smatrix2);
+
+ protected:
+ VectorView m_data;
+
+}; //end ConstBaseSymmetricMatrix
+
+
+template <size_t N>
+class BaseSymmetricMatrix : public ConstBaseSymmetricMatrix<N>
+{
+ public:
+ typedef ConstBaseSymmetricMatrix<N> base_type;
+
+
+ public:
+
+ BaseSymmetricMatrix (VectorView const& _data)
+ : base_type(_data)
+ {
+ }
+
+ double operator() (size_t i, size_t j) const
+ {
+ return m_data[base_type::get_index(i,j)];
+ }
+
+ double& operator() (size_t i, size_t j)
+ {
+ return m_data[base_type::get_index(i,j)];
+ }
+
+ template <size_t I, size_t J>
+ double& get()
+ {
+ BOOST_STATIC_ASSERT ((I < N && J < N));
+ return m_data[detail::index<I, J>::K];
+ }
+
+ void set_all (double x)
+ {
+ m_data.set_all(x);
+ }
+
+ SymmetricMatrixView<N-1> main_minor_view();
+
+ BaseSymmetricMatrix& transpose() const
+ {
+ return (*this);
+ }
+
+ BaseSymmetricMatrix& translate (double c)
+ {
+ m_data.translate(c);
+ return (*this);
+ }
+
+ BaseSymmetricMatrix& scale (double c)
+ {
+ m_data.scale(c);
+ return (*this);
+ }
+
+ BaseSymmetricMatrix& operator+= (base_type const& _smatrix)
+ {
+ m_data += (static_cast<const BaseSymmetricMatrix &>(_smatrix).m_data);
+ return (*this);
+ }
+
+ BaseSymmetricMatrix& operator-= (base_type const& _smatrix)
+ {
+ m_data -= (static_cast<const BaseSymmetricMatrix &>(_smatrix).m_data);
+ return (*this);
+ }
+
+ using base_type::DIM;
+ using base_type::DATA_SIZE;
+ using base_type::m_data;
+ using base_type::operator-;
+ using base_type::operator*;
+
+}; //end BaseSymmetricMatrix
+
+
+template <size_t N>
+class SymmetricMatrix : public BaseSymmetricMatrix<N>
+{
+ public:
+ typedef BaseSymmetricMatrix<N> base_type;
+ typedef typename base_type::base_type base_base_type;
+
+ using base_type::DIM;
+ using base_type::DATA_SIZE;
+ using base_type::m_data;
+
+ public:
+ SymmetricMatrix ()
+ : base_type (VectorView(m_adata, DATA_SIZE))
+ {
+ }
+
+ explicit
+ SymmetricMatrix (ConstVectorView _data)
+ : base_type (VectorView(m_adata, DATA_SIZE))
+ {
+ assert (_data.size() == DATA_SIZE);
+ m_data = _data;
+ }
+
+ explicit
+ SymmetricMatrix (const double _data[DATA_SIZE])
+ : base_type (VectorView(m_adata, DATA_SIZE))
+ {
+ std::copy (_data, _data + DATA_SIZE, m_adata);
+ }
+
+ SymmetricMatrix (SymmetricMatrix const& _smatrix)
+ : base_type (VectorView(m_adata, DATA_SIZE))
+ {
+ m_data = _smatrix.m_data;
+ }
+
+ explicit
+ SymmetricMatrix (base_base_type const& _smatrix)
+ : base_type (VectorView(m_adata, DATA_SIZE))
+ {
+ m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data;
+ }
+
+ explicit
+ SymmetricMatrix (ConstMatrixView const& _matrix)
+ : base_type (VectorView(m_adata, DATA_SIZE))
+ {
+ assert (_matrix.rows() == N && _matrix.columns() == N);
+ for (size_t i = 0; i < N; ++i)
+ for (size_t j = 0; j <= i ; ++j)
+ (*this)(i,j) = _matrix(i,j);
+ }
+
+ SymmetricMatrix& operator= (SymmetricMatrix const& _smatrix)
+ {
+ m_data = _smatrix.m_data;
+ return (*this);
+ }
+
+ SymmetricMatrix& operator= (base_base_type const& _smatrix)
+ {
+
+ m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data;
+ return (*this);
+ }
+
+ SymmetricMatrix& operator= (ConstMatrixView const& _matrix)
+ {
+ assert (_matrix.rows() == N && _matrix.columns() == N);
+ for (size_t i = 0; i < N; ++i)
+ for (size_t j = 0; j <= i ; ++j)
+ (*this)(i,j) = _matrix(i,j);
+
+ return (*this);
+ }
+
+ // needed for accessing m_adata
+ friend class ConstSymmetricMatrixView<DIM>;
+ friend class SymmetricMatrixView<DIM>;
+ private:
+ double m_adata[DATA_SIZE];
+}; //end SymmetricMatrix
+
+
+template <size_t N>
+class ConstSymmetricMatrixView : public ConstBaseSymmetricMatrix<N>
+{
+ public:
+ typedef ConstBaseSymmetricMatrix<N> base_type;
+
+ using base_type::DIM;
+ using base_type::DATA_SIZE;
+ using base_type::m_data;
+
+
+ public:
+
+ explicit
+ ConstSymmetricMatrixView (ConstVectorView _data)
+ : base_type (const_vector_view_cast(_data))
+ {
+ assert (_data.size() == DATA_SIZE);
+ }
+
+ explicit
+ ConstSymmetricMatrixView (const double _data[DATA_SIZE])
+ : base_type (const_vector_view_cast (ConstVectorView (_data, DATA_SIZE)))
+ {
+ }
+
+ ConstSymmetricMatrixView (const ConstSymmetricMatrixView & _smatrix)
+ : base_type (_smatrix.m_data)
+ {
+ }
+
+ ConstSymmetricMatrixView (const base_type & _smatrix)
+ : base_type (static_cast<const ConstSymmetricMatrixView &>(_smatrix).m_data)
+ {
+ }
+
+}; //end SymmetricMatrix
+
+
+// declaration needed for friend clause
+template <size_t N>
+void swap_view(SymmetricMatrixView<N> & m1, SymmetricMatrixView<N> & m2);
+
+
+template <size_t N>
+class SymmetricMatrixView : public BaseSymmetricMatrix<N>
+{
+ public:
+ typedef BaseSymmetricMatrix<N> base_type;
+ typedef typename base_type::base_type base_base_type;
+
+ using base_type::DIM;
+ using base_type::DATA_SIZE;
+ using base_type::m_data;
+
+ public:
+
+ explicit
+ SymmetricMatrixView (VectorView _data)
+ : base_type (_data)
+ {
+ assert (_data.size() == DATA_SIZE);
+ }
+
+ explicit
+ SymmetricMatrixView (double _data[DATA_SIZE])
+ : base_type (VectorView (_data, DATA_SIZE))
+ {
+ }
+
+ SymmetricMatrixView (const SymmetricMatrixView & _smatrix)
+ : base_type (_smatrix.m_data)
+ {
+ }
+
+ SymmetricMatrixView (SymmetricMatrix<DIM> & _smatrix)
+ : base_type (VectorView (_smatrix.m_adata, DATA_SIZE))
+ {
+ }
+
+ SymmetricMatrixView& operator= (const SymmetricMatrixView & _smatrix)
+ {
+ m_data = _smatrix.m_data;
+ return (*this);
+ }
+
+ SymmetricMatrixView& operator= (const base_base_type & _smatrix)
+ {
+ m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data;
+ return (*this);
+ }
+
+ friend
+ void swap_view<N>(SymmetricMatrixView & m1, SymmetricMatrixView & m2);
+
+}; //end SymmetricMatrix
+
+
+
+
+/*
+ * class ConstBaseSymmetricMatrix methods
+ */
+
+template <size_t N>
+inline
+std::string ConstBaseSymmetricMatrix<N>::str() const
+{
+ std::ostringstream oss;
+ oss << (*this);
+ return oss.str();
+}
+
+template <size_t N>
+inline
+ConstSymmetricMatrixView<N-1>
+ConstBaseSymmetricMatrix<N>::main_minor_const_view() const
+{
+ ConstVectorView data(m_data.get_gsl_vector()->data, DATA_SIZE - DIM);
+ ConstSymmetricMatrixView<N-1> mm(data);
+ return mm;
+}
+
+template <size_t N>
+inline
+SymmetricMatrix<N> ConstBaseSymmetricMatrix<N>::operator- () const
+{
+ SymmetricMatrix<N> result;
+ for (size_t i = 0; i < DATA_SIZE; ++i)
+ {
+ result.m_data[i] = -m_data[i];
+ }
+ return result;
+}
+
+
+/*
+ * class ConstBaseSymmetricMatrix friend free functions
+ */
+
+template <size_t N>
+inline
+bool operator== (ConstBaseSymmetricMatrix<N> const& _smatrix1,
+ ConstBaseSymmetricMatrix<N> const& _smatrix2)
+{
+ return (_smatrix1.m_data == _smatrix2.m_data);
+}
+
+/*
+ * class ConstBaseSymmetricMatrix related free functions
+ */
+
+template< size_t N, class charT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os,
+ const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+ os << "[[" << _matrix(0,0);
+ for (size_t j = 1; j < N; ++j)
+ {
+ os << ", " << _matrix(0,j);
+ }
+ os << "]";
+ for (size_t i = 1; i < N; ++i)
+ {
+ os << "\n [" << _matrix(i,0);
+ for (size_t j = 1; j < N; ++j)
+ {
+ os << ", " << _matrix(i,j);
+ }
+ os << "]";
+ }
+ os << "]";
+ return os;
+}
+
+
+/*
+ * class ConstBaseSymmetricMatrix specialized methods
+ */
+
+template<>
+inline
+size_t ConstBaseSymmetricMatrix<2>::get_index (size_t i, size_t j)
+{
+ return (i+j);
+}
+
+template<>
+inline
+size_t ConstBaseSymmetricMatrix<3>::get_index (size_t i, size_t j)
+{
+ size_t k = i + j;
+ if (i == 2 || j == 2) ++k;
+ return k;
+}
+
+
+/*
+ * class BaseSymmetricMatrix methods
+ */
+
+template <size_t N>
+inline
+SymmetricMatrixView<N-1> BaseSymmetricMatrix<N>::main_minor_view()
+{
+ VectorView data(m_data.get_gsl_vector()->data, DATA_SIZE - DIM);
+ SymmetricMatrixView<N-1> mm(data);
+ return mm;
+}
+
+
+/*
+ * class SymmetricMatrixView friend free functions
+ */
+
+template <size_t N>
+inline
+void swap_view(SymmetricMatrixView<N> & m1, SymmetricMatrixView<N> & m2)
+{
+ swap_view(m1.m_data, m2.m_data);
+}
+
+} /* end namespace NL*/ } /* end namespace Geom*/
+
+
+
+
+#endif // _NL_SYMMETRIC_MATRIX_FS_H_
+
+
+
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/vector.h b/src/2geom/numeric/vector.h
new file mode 100644
index 0000000..f28289f
--- /dev/null
+++ b/src/2geom/numeric/vector.h
@@ -0,0 +1,594 @@
+/*
+ * 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
+ {
+ for ( size_t i = 0; i < size(); ++i )
+ {
+ if ( (*this)[i] <= 0 ) return false;
+ }
+ return true;
+ }
+
+ bool is_negative() const
+ {
+ for ( size_t i = 0; i < size(); ++i )
+ {
+ if ( (*this)[i] >= 0 ) return false;
+ }
+ return true;
+ }
+
+ 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());
+ using std::swap;
+ swap(v1.m_vector, v2.m_vector);
+}
+
+inline
+void swap_any(Vector & v1, Vector & v2)
+{
+ using std::swap;
+ swap(v1.m_vector, v2.m_vector);
+ 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() );
+ using std::swap;
+ swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too
+}
+
+inline
+const VectorView & const_vector_view_cast (const ConstVectorView & view)
+{
+ const detail::BaseVectorImpl & bvi
+ = static_cast<const detail::BaseVectorImpl &>(view);
+ const VectorView & vv = reinterpret_cast<const VectorView &>(bvi);
+ return vv;
+}
+
+inline
+VectorView & const_vector_view_cast (ConstVectorView & view)
+{
+ detail::BaseVectorImpl & bvi
+ = static_cast<detail::BaseVectorImpl &>(view);
+ VectorView & vv = reinterpret_cast<VectorView &>(bvi);
+ return vv;
+}
+
+
+} } // 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:fileencoding=utf-8:textwidth=99 :