From 97a4fc5999be46ca00af89472361235e854a6d12 Mon Sep 17 00:00:00 2001 From: "Jon A. Cruz" Date: Fri, 4 Jul 2008 06:38:30 +0000 Subject: EOL fixup (bzr r6149) --- src/2geom/numeric/fitting-model.h | 846 +++++++++++++-------------- src/2geom/numeric/fitting-tool.h | 1064 +++++++++++++++++----------------- src/2geom/numeric/linear_system.h | 276 ++++----- src/2geom/numeric/matrix.h | 1110 ++++++++++++++++++------------------ src/2geom/numeric/vector.h | 1130 ++++++++++++++++++------------------- 5 files changed, 2213 insertions(+), 2213 deletions(-) (limited to 'src/2geom/numeric') diff --git a/src/2geom/numeric/fitting-model.h b/src/2geom/numeric/fitting-model.h index 145be40e4..cc3113372 100644 --- a/src/2geom/numeric/fitting-model.h +++ b/src/2geom/numeric/fitting-model.h @@ -1,423 +1,423 @@ -/* - * Fitting Models for Geom Types - * - * Authors: - * Marco Cecchetti - * - * 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 : +/* + * 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 index edacc663a..d589d86e3 100644 --- a/src/2geom/numeric/fitting-tool.h +++ b/src/2geom/numeric/fitting-tool.h @@ -1,532 +1,532 @@ -/* - * Fitting Tools - * - * Authors: - * Marco Cecchetti - * - * 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 : +/* + * 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 5b516c9e6..dc2a1d7e0 100644 --- a/src/2geom/numeric/linear_system.h +++ b/src/2geom/numeric/linear_system.h @@ -1,138 +1,138 @@ -/* - * LinearSystem class wraps some gsl routines for solving linear systems - * - * Authors: - * Marco Cecchetti - * - * 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 : +/* + * 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.h b/src/2geom/numeric/matrix.h index 64557a6f1..156b6e9a2 100644 --- a/src/2geom/numeric/matrix.h +++ b/src/2geom/numeric/matrix.h @@ -1,555 +1,555 @@ -/* - * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines; - * "views" mimic the semantic of C++ references: any operation performed - * on a "view" is actually performed on the "viewed object" - * - * Authors: - * Marco Cecchetti - * - * 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 : +/* + * 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 43a39a1ac..3e53405f4 100644 --- a/src/2geom/numeric/vector.h +++ b/src/2geom/numeric/vector.h @@ -1,565 +1,565 @@ -/* - * Vector, VectorView, ConstVectorView classes wrap the gsl vector routines; - * "views" mimic the semantic of C++ references: any operation performed - * on a "view" is actually performed on the "viewed object" - * - * Authors: - * Marco Cecchetti - * - * 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 : +/* + * 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