From e07de7ce18ff4cc440dc8947d61825122afea062 Mon Sep 17 00:00:00 2001 From: "Johan B. C. Engelen" Date: Thu, 3 Jul 2008 20:06:40 +0000 Subject: update 2geom (svn rev1433) (bzr r6134) --- src/2geom/numeric/fitting-model.h | 423 +++++++++++++++++++++ src/2geom/numeric/fitting-tool.h | 532 ++++++++++++++++++++++++++ src/2geom/numeric/linear_system.h | 227 +++++++----- src/2geom/numeric/matrix.cpp | 115 ++++++ src/2geom/numeric/matrix.h | 762 +++++++++++++++++++++++++++----------- src/2geom/numeric/vector.h | 750 ++++++++++++++++++++++++++++--------- 6 files changed, 2328 insertions(+), 481 deletions(-) create mode 100644 src/2geom/numeric/fitting-model.h create mode 100644 src/2geom/numeric/fitting-tool.h create mode 100644 src/2geom/numeric/matrix.cpp (limited to 'src/2geom/numeric') diff --git a/src/2geom/numeric/fitting-model.h b/src/2geom/numeric/fitting-model.h new file mode 100644 index 000000000..145be40e4 --- /dev/null +++ b/src/2geom/numeric/fitting-model.h @@ -0,0 +1,423 @@ +/* + * Fitting Models for Geom Types + * + * Authors: + * Marco Cecchetti + * + * 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 +{ + 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 +{ + public: + LFMPoly(size_t degree) + : LFMPowerBasis(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 +{ + 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 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 +{ + 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 +{ + 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 +{ + 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 objects +class LFMD2SBasis + : public LinearFittingModel< double, Point, D2 > +{ + 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 & 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 +{ + 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 m_bc; +}; + + +// this model generates Bezier curves +template< unsigned int N > +class LFMBezierCurve + : public LinearFittingModel< double, Point, BezierCurve > +{ + 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 & 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(bx, by); + } + + private: + LFMBezier mob; +}; + +} // end namespace NL +} // end namespace Geom + + +#endif // _NL_FITTING_MODEL_H_ + + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/numeric/fitting-tool.h b/src/2geom/numeric/fitting-tool.h new file mode 100644 index 000000000..edacc663a --- /dev/null +++ b/src/2geom/numeric/fitting-tool.h @@ -0,0 +1,532 @@ +/* + * Fitting Tools + * + * Authors: + * Marco Cecchetti + * + * 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 + + +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 + : public lsf_base +{ +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 base_type; + + using base_type::m_model; + using base_type::m_psdinv_matrix; + using base_type::total_samples; + +public: + lsf_solution( 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 + + +// 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 + : public lsf_base +{ +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 base_type; + + using base_type::m_model; + using base_type::m_psdinv_matrix; + using base_type::total_samples; + +public: + lsf_solution( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples), + m_solution(_model.size(), 2) + { + } + + solution_type& result(std::vector 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 const& old_sample_values, + std::vector 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 + + + + +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 + : public lsf_solution +{ + public: + typedef ModelT model_type; + typedef typename model_type::parameter_type parameter_type; + typedef typename model_type::value_type value_type; + typedef lsf_solution 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( 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 + + +// fitting tool for partially known models +template< typename ModelT> +class lsf_with_fixed_terms + : public lsf_solution +{ + public: + typedef ModelT model_type; + typedef typename model_type::parameter_type parameter_type; + typedef typename model_type::value_type value_type; + typedef lsf_solution 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( 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() + { + if (m_vector_view != NULL) + { + delete m_vector_view; + } + } + + protected: + Vector m_vector; + VectorView* m_vector_view; + +}; // end class lsf_with_fixed_terms + + +} // 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 + : public detail::lsf_with_fixed_terms +{ + public: + typedef ModelT model_type; + typedef detail::lsf_with_fixed_terms 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( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples) + { + } +}; // end class least_squeares_fitter + + +template< typename ModelT> +class least_squeares_fitter + : public detail::lsf_with_fixed_terms +{ + public: + typedef ModelT model_type; + typedef detail::lsf_with_fixed_terms 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( 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 + + +template< typename ModelT> +class least_squeares_fitter + : public detail::lsf_with_fixed_terms +{ + public: + typedef ModelT model_type; + typedef detail::lsf_with_fixed_terms 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( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples) + { + } + + solution_type& result(std::vector 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 + + +} // end namespace NL +} // end namespace Geom + + + +#endif // _NL_FITTING_TOOL_H_ + + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/numeric/linear_system.h b/src/2geom/numeric/linear_system.h index 7699c5224..5b516c9e6 100644 --- a/src/2geom/numeric/linear_system.h +++ b/src/2geom/numeric/linear_system.h @@ -1,89 +1,138 @@ -#ifndef _NL_LINEAR_SYSTEM_H_ -#define _NL_LINEAR_SYSTEM_H_ - - -#include - -#include - -#include "2geom/numeric/matrix.h" -#include "2geom/numeric/vector.h" - - -namespace Geom { namespace NL { - - -class LinearSystem -{ -public: - LinearSystem(Matrix & _matrix, Vector & _vector) - : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns()) - { - } - - const Vector & LU_solve() - { - assert( matrix().rows() == matrix().columns() - && matrix().rows() == vector().size() ); - int s; - gsl_permutation * p = gsl_permutation_alloc(matrix().rows()); - gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s); - gsl_linalg_LU_solve( matrix().get_gsl_matrix(), - p, - vector().get_gsl_vector(), - m_solution.get_gsl_vector() - ); - gsl_permutation_free(p); - return solution(); - } - - const Vector & SV_solve() - { - assert( matrix().rows() >= matrix().columns() - && matrix().rows() == vector().size() ); - - gsl_matrix* U = matrix().get_gsl_matrix(); - gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns()); - gsl_vector* S = gsl_vector_alloc(matrix().columns()); - gsl_vector* work = gsl_vector_alloc(matrix().columns()); - - gsl_linalg_SV_decomp( U, V, S, work ); - - gsl_vector* b = vector().get_gsl_vector(); - gsl_vector* x = m_solution.get_gsl_vector(); - - gsl_linalg_SV_solve( U, V, S, b, x); - - gsl_matrix_free(V); - gsl_vector_free(S); - gsl_vector_free(work); - - return solution(); - } - - Matrix & matrix() - { - return m_matrix; - } - - Vector & vector() - { - return m_vector; - } - - const Vector & solution() const - { - return m_solution; - } - -private: - Matrix & m_matrix; - Vector & m_vector; - Vector m_solution; -}; - - -} } // end namespaces - - -#endif /*_NL_LINEAR_SYSTEM_H_*/ +/* + * LinearSystem class wraps some gsl routines for solving linear systems + * + * Authors: + * Marco Cecchetti + * + * 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 + +#include + +#include <2geom/numeric/matrix.h> +#include <2geom/numeric/vector.h> + + +namespace Geom { namespace NL { + + +class LinearSystem +{ +public: + LinearSystem(MatrixView & _matrix, VectorView & _vector) + : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns()) + { + } + + LinearSystem(Matrix & _matrix, Vector & _vector) + : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns()) + { + } + + const Vector & LU_solve() + { + assert( matrix().rows() == matrix().columns() + && matrix().rows() == vector().size() ); + int s; + gsl_permutation * p = gsl_permutation_alloc(matrix().rows()); + gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s); + gsl_linalg_LU_solve( matrix().get_gsl_matrix(), + p, + vector().get_gsl_vector(), + m_solution.get_gsl_vector() + ); + gsl_permutation_free(p); + return solution(); + } + + const Vector & SV_solve() + { + assert( matrix().rows() >= matrix().columns() + && matrix().rows() == vector().size() ); + + gsl_matrix* U = matrix().get_gsl_matrix(); + gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns()); + gsl_vector* S = gsl_vector_alloc(matrix().columns()); + gsl_vector* work = gsl_vector_alloc(matrix().columns()); + + gsl_linalg_SV_decomp( U, V, S, work ); + + gsl_vector* b = vector().get_gsl_vector(); + gsl_vector* x = m_solution.get_gsl_vector(); + + gsl_linalg_SV_solve( U, V, S, b, x); + + gsl_matrix_free(V); + gsl_vector_free(S); + gsl_vector_free(work); + + return solution(); + } + + MatrixView & matrix() + { + return m_matrix; + } + + VectorView & vector() + { + return m_vector; + } + + const Vector & solution() const + { + return m_solution; + } + +private: + MatrixView m_matrix; + VectorView m_vector; + Vector m_solution; +}; + + +} } // end namespaces + + +#endif /*_NL_LINEAR_SYSTEM_H_*/ + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/numeric/matrix.cpp b/src/2geom/numeric/matrix.cpp new file mode 100644 index 000000000..dde80dfb6 --- /dev/null +++ b/src/2geom/numeric/matrix.cpp @@ -0,0 +1,115 @@ +/* + * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines; + * "views" mimic the semantic of C++ references: any operation performed + * on a "view" is actually performed on the "viewed object" + * + * Authors: + * Marco Cecchetti + * + * 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 +#include + + + +namespace Geom { namespace NL { + + +Vector operator*( detail::BaseMatrixImpl const& A, + detail::BaseVectorImpl const& v ) +{ + assert(A.columns() == v.size()); + + Vector result(A.rows(), 0.0); + for (size_t i = 0; i < A.rows(); ++i) + for (size_t j = 0; j < A.columns(); ++j) + result[i] += A(i,j) * v[j]; + + return result; +} + +Matrix operator*( detail::BaseMatrixImpl const& A, + detail::BaseMatrixImpl const& B ) +{ + assert(A.columns() == B.rows()); + + Matrix C(A.rows(), B.columns(), 0.0); + for (size_t i = 0; i < C.rows(); ++i) + for (size_t j = 0; j < C.columns(); ++j) + for (size_t k = 0; k < A.columns(); ++k) + C(i,j) += A(i,k) * B(k, j); + + return C; +} + +Matrix pseudo_inverse(detail::BaseMatrixImpl const& A) +{ + + Matrix U(A); + Matrix V(A.columns(), A.columns()); + Vector s(A.columns()); + gsl_vector* work = gsl_vector_alloc(A.columns()); + + gsl_linalg_SV_decomp( U.get_gsl_matrix(), + V.get_gsl_matrix(), + s.get_gsl_vector(), + work ); + + Matrix P(A.columns(), A.rows(), 0.0); + + int sz = s.size(); + while ( sz-- > 0 && s[sz] == 0 ) {} + ++sz; + if (sz == 0) return P; + VectorView sv(s, sz); + + for (size_t i = 0; i < sv.size(); ++i) + { + VectorView v = V.column_view(i); + v.scale(1/sv[i]); + for (size_t h = 0; h < P.rows(); ++h) + for (size_t k = 0; k < P.columns(); ++k) + P(h,k) += V(h,i) * U(k,i); + } + + return P; +} + +} } // end namespaces + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/numeric/matrix.h b/src/2geom/numeric/matrix.h index bdd647e53..64557a6f1 100644 --- a/src/2geom/numeric/matrix.h +++ b/src/2geom/numeric/matrix.h @@ -1,207 +1,555 @@ -#ifndef _NL_MATRIX_H_ -#define _NL_MATRIX_H_ - - -#include -#include -#include - -#include - - - -namespace Geom { namespace NL { - -class Matrix; -void swap(Matrix & m1, Matrix & m2); - - -class Matrix -{ -public: - // the matrix is not inizialized - Matrix(size_t n1, size_t n2) - : m_rows(n1), m_columns(n2) - { - m_matrix = gsl_matrix_alloc(n1, n2); - } - - Matrix(size_t n1, size_t n2, double x) - :m_rows(n1), m_columns(n2) - { - m_matrix = gsl_matrix_alloc(n1, n2); - gsl_matrix_set_all(m_matrix, x); - } - - Matrix( Matrix const& _matrix ) - : m_rows(_matrix.rows()), m_columns(_matrix.columns()) - { - m_matrix = gsl_matrix_alloc(rows(), columns()); - gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); - } - -// explicit -// Matrix( gsl_matrix* m, size_t n1, size_t n2) -// : m_rows(n1), m_columns(n2), m_matrix(m) -// { -// } - - Matrix & operator=(Matrix const& _matrix) - { - assert( rows() == _matrix.rows() && columns() == _matrix.columns() ); - gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); - - return *this; - } - - virtual ~Matrix() - { - gsl_matrix_free(m_matrix); - } - - void set_all( double x = 0) - { - gsl_matrix_set_all(m_matrix, x); - } - - void set_identity() - { - gsl_matrix_set_identity(m_matrix); - } - - const double & operator() (size_t i, size_t j) const - { - return *gsl_matrix_const_ptr(m_matrix, i, j); - } - - double & operator() (size_t i, size_t j) - { - return *gsl_matrix_ptr(m_matrix, i, j); - } - - gsl_matrix* get_gsl_matrix() - { - return m_matrix; - } - - void swap_rows(size_t i, size_t j) - { - gsl_matrix_swap_rows(m_matrix, i, j); - } - - void swap_columns(size_t i, size_t j) - { - gsl_matrix_swap_columns(m_matrix, i, j); - } - - Matrix & transpose() - { - gsl_matrix_transpose(m_matrix); - return (*this); - } - - Matrix & scale(double x) - { - gsl_matrix_scale(m_matrix, x); - return (*this); - } - - Matrix & translate(double x) - { - gsl_matrix_add_constant(m_matrix, x); - return (*this); - } - - Matrix & operator+=(Matrix const& _matrix) - { - gsl_matrix_add(m_matrix, _matrix.m_matrix); - return (*this); - } - - Matrix & operator-=(Matrix const& _matrix) - { - gsl_matrix_sub(m_matrix, _matrix.m_matrix); - return (*this); - } - - bool is_zero() const - { - return gsl_matrix_isnull(m_matrix); - } - - bool is_positive() const - { - return gsl_matrix_ispos(m_matrix); - } - - bool is_negative() const - { - return gsl_matrix_isneg(m_matrix); - } - - bool is_non_negative() const - { - for ( unsigned int i = 0; i < rows(); ++i ) - { - for ( unsigned int j = 0; j < columns(); ++j ) - { - if ( (*this)(i,j) < 0 ) return false; - } - } - return true; - } - - double max() const - { - return gsl_matrix_max(m_matrix); - } - - double min() const - { - return gsl_matrix_min(m_matrix); - } - - std::pair - max_index() const - { - std::pair indices; - gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second)); - return indices; - } - - std::pair - min_index() const - { - std::pair indices; - gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second)); - return indices; - } - - friend - void swap(Matrix & m1, Matrix & m2); - - size_t rows() const - { - return m_rows; - } - - size_t columns() const - { - return m_columns; - } - -private: - size_t m_rows, m_columns; - gsl_matrix* m_matrix; -}; - -void swap(Matrix & m1, Matrix & m2) -{ - assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() ); - std::swap(m1.m_matrix, m2.m_matrix); -} - - -} } // end namespaces - -#endif /*_NL_MATRIX_H_*/ +/* + * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines; + * "views" mimic the semantic of C++ references: any operation performed + * on a "view" is actually performed on the "viewed object" + * + * Authors: + * Marco Cecchetti + * + * 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 +#include // for std::pair +#include // for std::swap +#include +#include + +#include +#include + + +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 + max_index() const + { + std::pair indices; + gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second)); + return indices; + } + + std::pair + min_index() const + { + std::pair 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 & +operator<< (std::basic_ostream & 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( base_type::transpose() ); + } + + Matrix & scale(double x) + { + return static_cast( base_type::scale(x) ); + } + + Matrix & translate(double x) + { + return static_cast( base_type::translate(x) ); + } + + Matrix & operator+=(base_type::base_type const& _matrix) + { + return static_cast( base_type::operator+=(_matrix) ); + } + + Matrix & operator-=(base_type::base_type const& _matrix) + { + return static_cast( 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( &(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( &(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( &(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( base_type::transpose() ); + } + + MatrixView & scale(double x) + { + return static_cast( base_type::scale(x) ); + } + + MatrixView & translate(double x) + { + return static_cast( base_type::translate(x) ); + } + + MatrixView & operator+=(base_type::base_type const& _matrix) + { + return static_cast( base_type::operator+=(_matrix) ); + } + + MatrixView & operator-=(base_type::base_type const& _matrix) + { + return static_cast( base_type::operator-=(_matrix) ); + } + + friend + void swap_view(MatrixView & m1, MatrixView & m2); + + private: + gsl_matrix_view m_matrix_view; + +}; // end class MatrixView + + +inline +void swap_view(MatrixView & m1, MatrixView & m2) +{ + assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() ); + std::swap(m1.m_matrix_view, m2.m_matrix_view); +} + +Vector operator*( detail::BaseMatrixImpl const& A, + detail::BaseVectorImpl const& v ); + +Matrix operator*( detail::BaseMatrixImpl const& A, + detail::BaseMatrixImpl const& B ); + +Matrix pseudo_inverse(detail::BaseMatrixImpl const& A); + +} } // end namespaces + +#endif /*_NL_MATRIX_H_*/ + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/numeric/vector.h b/src/2geom/numeric/vector.h index b25861e22..43a39a1ac 100644 --- a/src/2geom/numeric/vector.h +++ b/src/2geom/numeric/vector.h @@ -1,185 +1,565 @@ -#ifndef _NL_VECTOR_H_ -#define _NL_VECTOR_H_ - -#include -#include - -#include - - -namespace Geom { namespace NL { - -class Vector; -void swap(Vector & v1, Vector & v2); - - -class Vector -{ -public: - Vector(size_t n) - : m_size(n) - { - m_vector = gsl_vector_alloc(n); - } - - Vector(size_t n, double x) - : m_size(n) - { - m_vector = gsl_vector_alloc(n); - gsl_vector_set_all(m_vector, x); - } - - // create a vector with n elements all set to zero - // but the i-th that is set to 1 - Vector(size_t n, size_t i) - : m_size(n) - { - m_vector = gsl_vector_alloc(n); - gsl_vector_set_basis(m_vector, i); - } - - Vector(Vector const& _vector) - : m_size(_vector.size()) - { - m_vector = gsl_vector_alloc(size()); - gsl_vector_memcpy(m_vector, _vector.m_vector); - } - - virtual ~Vector() - { - gsl_vector_free(m_vector); - } - - Vector & operator=(Vector const& _vector) - { - assert( size() == _vector.size() ); - gsl_vector_memcpy(m_vector, _vector.m_vector); - return (*this); - } - - void set_all(double x = 0) - { - gsl_vector_set_all(m_vector, x); - } - - void set_basis(size_t i) - { - gsl_vector_set_basis(m_vector, i); - } - - double const& operator[](size_t i) const - { - return *gsl_vector_const_ptr(m_vector, i); - } - - double & operator[](size_t i) - { - return *gsl_vector_ptr(m_vector, i); - } - - gsl_vector* get_gsl_vector() - { - return m_vector; - } - - void swap_elements(size_t i, size_t j) - { - gsl_vector_swap_elements(m_vector, i, j); - } - - void reverse() - { - gsl_vector_reverse(m_vector); - } - - Vector & scale(double x) - { - gsl_vector_scale(m_vector, x); - return (*this); - } - - Vector & translate(double x) - { - gsl_vector_add_constant(m_vector, x); - return (*this); - } - - Vector & operator+=(Vector const& _vector) - { - gsl_vector_add(m_vector, _vector.m_vector); - return (*this); - } - - Vector & operator-=(Vector const& _vector) - { - gsl_vector_sub(m_vector, _vector.m_vector); - return (*this); - } - - bool is_zero() const - { - return gsl_vector_isnull(m_vector); - } - - bool is_positive() const - { - return gsl_vector_ispos(m_vector); - } - - bool is_negative() const - { - return gsl_vector_isneg(m_vector); - } - - bool is_non_negative() const - { - for ( size_t i = 0; i < size(); ++i ) - { - if ( (*this)[i] < 0 ) return false; - } - return true; - } - - double max() const - { - return gsl_vector_max(m_vector); - } - - double min() const - { - return gsl_vector_min(m_vector); - } - - size_t max_index() const - { - return gsl_vector_max_index(m_vector); - } - - size_t min_index() const - { - return gsl_vector_min_index(m_vector); - } - - friend - void swap(Vector & v1, Vector & v2); - - size_t size() const - { - return m_size; - } - -private: - size_t m_size; - gsl_vector* m_vector; -}; - -void swap(Vector & v1, Vector & v2) -{ - assert( v1.size() == v2.size() ); - std::swap(v1.m_vector, v2.m_vector); -} - -} } // end namespaces - - -#endif /*_NL_VECTOR_H_*/ +/* + * Vector, VectorView, ConstVectorView classes wrap the gsl vector routines; + * "views" mimic the semantic of C++ references: any operation performed + * on a "view" is actually performed on the "viewed object" + * + * Authors: + * Marco Cecchetti + * + * 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 +#include // for std::swap +#include +#include +#include + + +#include +#include + + +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 & +operator<< (std::basic_ostream & 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( base_type::scale(x) ); + } + + Vector & translate(double x) + { + return static_cast( base_type::translate(x) ); + } + + Vector & operator+=(base_type::base_type const& _vector) + { + return static_cast( base_type::operator+=(_vector) ); + } + + Vector & operator-=(base_type::base_type const& _vector) + { + return static_cast( 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( &(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( &(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( &(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( &(m_vector_view.vector) ); + } + + explicit + ConstVectorView(gsl_vector_const_view _gsl_vector_view) + : m_vector_view(_gsl_vector_view) + { + m_vector = const_cast( &(m_vector_view.vector) ); + m_size = m_vector->size; + } + + explicit + ConstVectorView(const std::vector& _vector) + : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) ) + { + m_vector = const_cast( &(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( &(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( &(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 & _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( base_type::scale(x) ); + } + + VectorView & translate(double x) + { + return static_cast( base_type::translate(x) ); + } + + VectorView & operator+=(base_type::base_type const& _vector) + { + return static_cast( base_type::operator+=(_vector) ); + } + + VectorView & operator-=(base_type::base_type const& _vector) + { + return static_cast( 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 : -- cgit v1.2.3