diff options
Diffstat (limited to 'include/2geom/numeric')
| -rw-r--r-- | include/2geom/numeric/fitting-model.h | 521 | ||||
| -rw-r--r-- | include/2geom/numeric/fitting-tool.h | 562 | ||||
| -rw-r--r-- | include/2geom/numeric/linear_system.h | 138 | ||||
| -rw-r--r-- | include/2geom/numeric/matrix.h | 603 | ||||
| -rw-r--r-- | include/2geom/numeric/symmetric-matrix-fs-operation.h | 102 | ||||
| -rw-r--r-- | include/2geom/numeric/symmetric-matrix-fs-trace.h | 427 | ||||
| -rw-r--r-- | include/2geom/numeric/symmetric-matrix-fs.h | 733 | ||||
| -rw-r--r-- | include/2geom/numeric/vector.h | 594 |
8 files changed, 3680 insertions, 0 deletions
diff --git a/include/2geom/numeric/fitting-model.h b/include/2geom/numeric/fitting-model.h new file mode 100644 index 0000000..0316f57 --- /dev/null +++ b/include/2geom/numeric/fitting-model.h @@ -0,0 +1,521 @@ +/* + * Fitting Models for Geom Types + * + * Authors: + * Marco Cecchetti <mrcekets at gmail.com> + * + * Copyright 2008 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + + +#ifndef _NL_FITTING_MODEL_H_ +#define _NL_FITTING_MODEL_H_ + + +#include <2geom/d2.h> +#include <2geom/sbasis.h> +#include <2geom/bezier.h> +#include <2geom/bezier-curve.h> +#include <2geom/polynomial.h> +#include <2geom/ellipse.h> +#include <2geom/circle.h> +#include <2geom/utils.h> +#include <2geom/conicsec.h> + + +namespace Geom { namespace NL { + +/* + * A model is an abstraction for an expression dependent from a parameter where + * the coefficients of this expression are the unknowns of the fitting problem. + * For a ceratain number of parameter values we know the related values + * the expression evaluates to: from each parameter value we get a row of + * the matrix of the fitting problem, from each expression value we get + * the related constant term. + * Example: given the model a*x^2 + b*x + c = 0; from x = 1 we get + * the equation a + b + c = 0, in this example the constant term is always + * the same for each parameter value. + * + * A model is required to implement 3 methods: + * + * - size : returns the number of unknown coefficients that appear in + * the expression of the fitting problem; + * - feed : its input is a parameter value and the related expression value, + * it generates a matrix row and a new entry of the constant vector + * of the fitting problem; + * - instance : it has an input parameter represented by the raw vector + * solution of the fitting problem and an output parameter + * of type InstanceType that return a specific object that is + * generated using the fitting problem solution, in the example + * above the object could be a Poly type. + */ + +/* + * completely unknown models must inherit from this template class; + * example: the model a*x^2 + b*x + c = 0 to be solved wrt a, b, c; + * example: the model A(t) = known_sample_value_at(t) to be solved wrt + * the coefficients of the curve A(t) expressed in S-Basis form; + * parameter type: the type of x and t variable in the examples above; + * value type: the type of the known sample values (in the first example + * is constant ) + * instance type: the type of the objects produced by using + * the fitting raw data solution + */ + + + + +template< typename ParameterType, typename ValueType, typename InstanceType > +class LinearFittingModel +{ + public: + typedef ParameterType parameter_type; + typedef ValueType value_type; + typedef InstanceType instance_type; + + static const bool WITH_FIXED_TERMS = false; + + /* + * a LinearFittingModel must implement the following methods: + * + * void feed( VectorView & vector, + * parameter_type const& sample_parameter ) const; + * + * size_t size() const; + * + * void instance(instance_type &, raw_type const& raw_data) const; + * + */ +}; + + +/* + * partially known models must inherit from this template class + * example: the model a*x^2 + 2*x + c = 0 to be solved wrt a and c + */ +template< typename ParameterType, typename ValueType, typename InstanceType > +class LinearFittingModelWithFixedTerms +{ + public: + typedef ParameterType parameter_type; + typedef ValueType value_type; + typedef InstanceType instance_type; + + static const bool WITH_FIXED_TERMS = true; + + /* + * a LinearFittingModelWithFixedTerms must implement the following methods: + * + * void feed( VectorView & vector, + * value_type & fixed_term, + * parameter_type const& sample_parameter ) const; + * + * size_t size() const; + * + * void instance(instance_type &, raw_type const& raw_data) const; + * + */ + + +}; + + +// incomplete model, it can be inherited to make up different kinds of +// instance type; the raw data is a vector of coefficients of a polynomial +// represented in standard power basis +template< typename InstanceType > +class LFMPowerBasis + : public LinearFittingModel<double, double, InstanceType> +{ + public: + LFMPowerBasis(size_t degree) + : m_size(degree + 1) + { + } + + void feed( VectorView & coeff, double sample_parameter ) const + { + coeff[0] = 1; + double x_i = 1; + for (size_t i = 1; i < coeff.size(); ++i) + { + x_i *= sample_parameter; + coeff[i] = x_i; + } + } + + size_t size() const + { + return m_size; + } + + private: + size_t m_size; +}; + + +// this model generates Geom::Poly objects +class LFMPoly + : public LFMPowerBasis<Poly> +{ + public: + LFMPoly(size_t degree) + : LFMPowerBasis<Poly>(degree) + { + } + + void instance(Poly & poly, ConstVectorView const& raw_data) const + { + poly.clear(); + poly.resize(size()); + for (size_t i = 0; i < raw_data.size(); ++i) + { + poly[i] = raw_data[i]; + } + } +}; + + +// incomplete model, it can be inherited to make up different kinds of +// instance type; the raw data is a vector of coefficients of a polynomial +// represented in standard power basis with leading term coefficient equal to 1 +template< typename InstanceType > +class LFMNormalizedPowerBasis + : public LinearFittingModelWithFixedTerms<double, double, InstanceType> +{ + public: + LFMNormalizedPowerBasis(size_t _degree) + : m_model( _degree - 1) + { + assert(_degree > 0); + } + + + void feed( VectorView & coeff, + double & known_term, + double sample_parameter ) const + { + m_model.feed(coeff, sample_parameter); + known_term = coeff[m_model.size()-1] * sample_parameter; + } + + size_t size() const + { + return m_model.size(); + } + + private: + LFMPowerBasis<InstanceType> m_model; +}; + + +// incomplete model, it can be inherited to make up different kinds of +// instance type; the raw data is a vector of coefficients of the equation +// of an ellipse curve +//template< typename InstanceType > +//class LFMEllipseEquation +// : public LinearFittingModelWithFixedTerms<Point, double, InstanceType> +//{ +// public: +// void feed( VectorView & coeff, double & fixed_term, Point const& p ) const +// { +// coeff[0] = p[X] * p[Y]; +// coeff[1] = p[Y] * p[Y]; +// coeff[2] = p[X]; +// coeff[3] = p[Y]; +// coeff[4] = 1; +// fixed_term = p[X] * p[X]; +// } +// +// size_t size() const +// { +// return 5; +// } +//}; + +// incomplete model, it can be inherited to make up different kinds of +// instance type; the raw data is a vector of coefficients of the equation +// of a conic section +template< typename InstanceType > +class LFMConicEquation + : public LinearFittingModelWithFixedTerms<Point, double, InstanceType> +{ + public: + void feed( VectorView & coeff, double & fixed_term, Point const& p ) const + { + coeff[0] = p[X] * p[Y]; + coeff[1] = p[Y] * p[Y]; + coeff[2] = p[X]; + coeff[3] = p[Y]; + coeff[4] = 1; + fixed_term = p[X] * p[X]; + } + + size_t size() const + { + return 5; + } +}; + +// this model generates Ellipse curves +class LFMConicSection + : public LFMConicEquation<xAx> +{ + public: + void instance(xAx & c, ConstVectorView const& coeff) const + { + c.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]); + } +}; + +// this model generates Ellipse curves +class LFMEllipse + : public LFMConicEquation<Ellipse> +{ + public: + void instance(Ellipse & e, ConstVectorView const& coeff) const + { + e.setCoefficients(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]); + } +}; + + +// incomplete model, it can be inherited to make up different kinds of +// instance type; the raw data is a vector of coefficients of the equation +// of a circle curve +template< typename InstanceType > +class LFMCircleEquation + : public LinearFittingModelWithFixedTerms<Point, double, InstanceType> +{ + public: + void feed( VectorView & coeff, double & fixed_term, Point const& p ) const + { + coeff[0] = p[X]; + coeff[1] = p[Y]; + coeff[2] = 1; + fixed_term = p[X] * p[X] + p[Y] * p[Y]; + } + + size_t size() const + { + return 3; + } +}; + + +// this model generates Ellipse curves +class LFMCircle + : public LFMCircleEquation<Circle> +{ + public: + void instance(Circle & c, ConstVectorView const& coeff) const + { + c.setCoefficients(1, coeff[0], coeff[1], coeff[2]); + } +}; + + +// this model generates SBasis objects +class LFMSBasis + : public LinearFittingModel<double, double, SBasis> +{ + public: + LFMSBasis( size_t _order ) + : m_size( 2*(_order+1) ), + m_order(_order) + { + } + + void feed( VectorView & coeff, double t ) const + { + double u0 = 1-t; + double u1 = t; + double s = u0 * u1; + coeff[0] = u0; + coeff[1] = u1; + for (size_t i = 2; i < size(); i+=2) + { + u0 *= s; + u1 *= s; + coeff[i] = u0; + coeff[i+1] = u1; + } + } + + size_t size() const + { + return m_size; + } + + void instance(SBasis & sb, ConstVectorView const& raw_data) const + { + sb.resize(m_order+1); + for (unsigned int i = 0, k = 0; i < raw_data.size(); i+=2, ++k) + { + sb[k][0] = raw_data[i]; + sb[k][1] = raw_data[i+1]; + } + } + + private: + size_t m_size; + size_t m_order; +}; + + +// this model generates D2<SBasis> objects +class LFMD2SBasis + : public LinearFittingModel< double, Point, D2<SBasis> > +{ + public: + LFMD2SBasis( size_t _order ) + : mosb(_order) + { + } + + void feed( VectorView & coeff, double t ) const + { + mosb.feed(coeff, t); + } + + size_t size() const + { + return mosb.size(); + } + + void instance(D2<SBasis> & d2sb, ConstMatrixView const& raw_data) const + { + mosb.instance(d2sb[X], raw_data.column_const_view(X)); + mosb.instance(d2sb[Y], raw_data.column_const_view(Y)); + } + + private: + LFMSBasis mosb; +}; + + +// this model generates Bezier objects +class LFMBezier + : public LinearFittingModel<double, double, Bezier> +{ + public: + LFMBezier( size_t _order ) + : m_size(_order + 1), + m_order(_order) + { + binomial_coefficients(m_bc, m_order); + } + + void feed( VectorView & coeff, double t ) const + { + double s = 1; + for (size_t i = 0; i < size(); ++i) + { + coeff[i] = s * m_bc[i]; + s *= t; + } + double u = 1-t; + s = 1; + for (size_t i = size()-1; i > 0; --i) + { + coeff[i] *= s; + s *= u; + } + coeff[0] *= s; + } + + size_t size() const + { + return m_size; + } + + void instance(Bezier & b, ConstVectorView const& raw_data) const + { + assert(b.size() == raw_data.size()); + for (unsigned int i = 0; i < raw_data.size(); ++i) + { + b[i] = raw_data[i]; + } + } + + private: + size_t m_size; + size_t m_order; + std::vector<size_t> m_bc; +}; + + +// this model generates Bezier curves +template <unsigned degree> +class LFMBezierCurveN + : public LinearFittingModel< double, Point, BezierCurveN<degree> > +{ + public: + LFMBezierCurveN() + : mob(degree+1) + { + } + + void feed( VectorView & coeff, double t ) const + { + mob.feed(coeff, t); + } + + size_t size() const + { + return mob.size(); + } + + void instance(BezierCurveN<degree> & bc, ConstMatrixView const& raw_data) const + { + Bezier bx(degree); + Bezier by(degree); + mob.instance(bx, raw_data.column_const_view(X)); + mob.instance(by, raw_data.column_const_view(Y)); + bc = BezierCurveN<degree>(bx, by); + } + + private: + LFMBezier mob; +}; + +} // end namespace NL +} // end namespace Geom + + +#endif // _NL_FITTING_MODEL_H_ + + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : diff --git a/include/2geom/numeric/fitting-tool.h b/include/2geom/numeric/fitting-tool.h new file mode 100644 index 0000000..78d66ca --- /dev/null +++ b/include/2geom/numeric/fitting-tool.h @@ -0,0 +1,562 @@ +/* + * Fitting Tools + * + * Authors: + * Marco Cecchetti <mrcekets at gmail.com> + * + * Copyright 2008 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + + +#ifndef _NL_FITTING_TOOL_H_ +#define _NL_FITTING_TOOL_H_ + + +#include <2geom/numeric/vector.h> +#include <2geom/numeric/matrix.h> + +#include <2geom/point.h> + +#include <vector> + + +/* + * The least_square_fitter class represents a tool for solving a fitting + * problem with respect to a given model that represents an expression + * dependent from a parameter where the coefficients of this expression + * are the unknowns of the fitting problem. + * The minimizing solution is found by computing the pseudo-inverse + * of the problem matrix + */ + + +namespace Geom { namespace NL { + +namespace detail { + +template< typename ModelT> +class lsf_base +{ + public: + typedef ModelT model_type; + typedef typename model_type::parameter_type parameter_type; + typedef typename model_type::value_type value_type; + + lsf_base( model_type const& _model, size_t forecasted_samples ) + : m_model(_model), + m_total_samples(0), + m_matrix(forecasted_samples, m_model.size()), + m_psdinv_matrix(NULL) + { + } + + // compute pseudo inverse + void update() + { + if (total_samples() == 0) return; + if (m_psdinv_matrix != NULL) + { + delete m_psdinv_matrix; + } + MatrixView mv(m_matrix, 0, 0, total_samples(), m_matrix.columns()); + m_psdinv_matrix = new Matrix( pseudo_inverse(mv) ); + assert(m_psdinv_matrix != NULL); + } + + size_t total_samples() const + { + return m_total_samples; + } + + bool is_full() const + { + return (total_samples() == m_matrix.rows()); + } + + void clear() + { + m_total_samples = 0; + } + + virtual + ~lsf_base() + { + if (m_psdinv_matrix != NULL) + { + delete m_psdinv_matrix; + } + } + + protected: + const model_type & m_model; + size_t m_total_samples; + Matrix m_matrix; + Matrix* m_psdinv_matrix; + +}; // end class lsf_base + + + + +template< typename ModelT, typename ValueType = typename ModelT::value_type> +class lsf_solution +{ +}; + +// a fitting process on samples with value of type double +// produces a solution of type Vector +template< typename ModelT> +class lsf_solution<ModelT, double> + : public lsf_base<ModelT> +{ +public: + typedef ModelT model_type; + typedef typename model_type::parameter_type parameter_type; + typedef typename model_type::value_type value_type; + typedef Vector solution_type; + typedef lsf_base<model_type> base_type; + + using base_type::m_model; + using base_type::m_psdinv_matrix; + using base_type::total_samples; + +public: + lsf_solution<ModelT, double>( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples), + m_solution(_model.size()) + { + } + + template< typename VectorT > + solution_type& result(VectorT const& sample_values) + { + assert(sample_values.size() == total_samples()); + ConstVectorView sv(sample_values); + m_solution = (*m_psdinv_matrix) * sv; + return m_solution; + } + + // a comparison between old sample values and the new ones is performed + // in order to minimize computation + // prerequisite: + // old_sample_values.size() == new_sample_values.size() + // no update() call can be performed between two result invocations + template< typename VectorT > + solution_type& result( VectorT const& old_sample_values, + VectorT const& new_sample_values ) + { + assert(old_sample_values.size() == total_samples()); + assert(new_sample_values.size() == total_samples()); + Vector diff(total_samples()); + for (size_t i = 0; i < diff.size(); ++i) + { + diff[i] = new_sample_values[i] - old_sample_values[i]; + } + Vector column(m_model.size()); + Vector delta(m_model.size(), 0.0); + for (size_t i = 0; i < diff.size(); ++i) + { + if (diff[i] != 0) + { + column = m_psdinv_matrix->column_view(i); + column.scale(diff[i]); + delta += column; + } + } + m_solution += delta; + return m_solution; + } + + solution_type& result() + { + return m_solution; + } + +private: + solution_type m_solution; + +}; // end class lsf_solution<ModelT, double> + + +// a fitting process on samples with value of type Point +// produces a solution of type Matrix (with 2 columns) +template< typename ModelT> +class lsf_solution<ModelT, Point> + : public lsf_base<ModelT> +{ +public: + typedef ModelT model_type; + typedef typename model_type::parameter_type parameter_type; + typedef typename model_type::value_type value_type; + typedef Matrix solution_type; + typedef lsf_base<model_type> base_type; + + using base_type::m_model; + using base_type::m_psdinv_matrix; + using base_type::total_samples; + +public: + lsf_solution<ModelT, Point>( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples), + m_solution(_model.size(), 2) + { + } + + solution_type& result(std::vector<Point> const& sample_values) + { + assert(sample_values.size() == total_samples()); + Matrix svm(total_samples(), 2); + for (size_t i = 0; i < total_samples(); ++i) + { + svm(i, X) = sample_values[i][X]; + svm(i, Y) = sample_values[i][Y]; + } + m_solution = (*m_psdinv_matrix) * svm; + return m_solution; + } + + // a comparison between old sample values and the new ones is performed + // in order to minimize computation + // prerequisite: + // old_sample_values.size() == new_sample_values.size() + // no update() call can to be performed between two result invocations + solution_type& result( std::vector<Point> const& old_sample_values, + std::vector<Point> const& new_sample_values ) + { + assert(old_sample_values.size() == total_samples()); + assert(new_sample_values.size() == total_samples()); + Matrix diff(total_samples(), 2); + for (size_t i = 0; i < total_samples(); ++i) + { + diff(i, X) = new_sample_values[i][X] - old_sample_values[i][X]; + diff(i, Y) = new_sample_values[i][Y] - old_sample_values[i][Y]; + } + Vector column(m_model.size()); + Matrix delta(m_model.size(), 2, 0.0); + VectorView deltax = delta.column_view(X); + VectorView deltay = delta.column_view(Y); + for (size_t i = 0; i < total_samples(); ++i) + { + if (diff(i, X) != 0) + { + column = m_psdinv_matrix->column_view(i); + column.scale(diff(i, X)); + deltax += column; + } + if (diff(i, Y) != 0) + { + column = m_psdinv_matrix->column_view(i); + column.scale(diff(i, Y)); + deltay += column; + } + } + m_solution += delta; + return m_solution; + } + + solution_type& result() + { + return m_solution; + } + +private: + solution_type m_solution; + +}; // end class lsf_solution<ModelT, Point> + + + + +template< typename ModelT, + bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS > +class lsf_with_fixed_terms +{ +}; + + +// fitting tool for completely unknown models +template< typename ModelT> +class lsf_with_fixed_terms<ModelT, false> + : public lsf_solution<ModelT> +{ + public: + typedef ModelT model_type; + typedef typename model_type::parameter_type parameter_type; + typedef typename model_type::value_type value_type; + typedef lsf_solution<model_type> base_type; + typedef typename base_type::solution_type solution_type; + + using base_type::total_samples; + using base_type::is_full; + using base_type::m_matrix; + using base_type::m_total_samples; + using base_type::m_model; + + public: + lsf_with_fixed_terms<ModelT, false>( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples) + { + } + + void append(parameter_type const& sample_parameter) + { + assert(!is_full()); + VectorView row = m_matrix.row_view(total_samples()); + m_model.feed(row, sample_parameter); + ++m_total_samples; + } + + void append_copy(size_t sample_index) + { + assert(!is_full()); + assert(sample_index < total_samples()); + VectorView dest_row = m_matrix.row_view(total_samples()); + VectorView source_row = m_matrix.row_view(sample_index); + dest_row = source_row; + ++m_total_samples; + } + +}; // end class lsf_with_fixed_terms<ModelT, false> + + +// fitting tool for partially known models +template< typename ModelT> +class lsf_with_fixed_terms<ModelT, true> + : public lsf_solution<ModelT> +{ + public: + typedef ModelT model_type; + typedef typename model_type::parameter_type parameter_type; + typedef typename model_type::value_type value_type; + typedef lsf_solution<model_type> base_type; + typedef typename base_type::solution_type solution_type; + + using base_type::total_samples; + using base_type::is_full; + using base_type::m_matrix; + using base_type::m_total_samples; + using base_type::m_model; + + public: + lsf_with_fixed_terms<ModelT, true>( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples), + m_vector(forecasted_samples), + m_vector_view(NULL) + { + } + void append(parameter_type const& sample_parameter) + { + assert(!is_full()); + VectorView row = m_matrix.row_view(total_samples()); + m_model.feed(row, m_vector[total_samples()], sample_parameter); + ++m_total_samples; + } + + void append_copy(size_t sample_index) + { + assert(!is_full()); + assert(sample_index < total_samples()); + VectorView dest_row = m_matrix.row_view(total_samples()); + VectorView source_row = m_matrix.row_view(sample_index); + dest_row = source_row; + m_vector[total_samples()] = m_vector[sample_index]; + ++m_total_samples; + } + + void update() + { + base_type::update(); + if (total_samples() == 0) return; + if (m_vector_view != NULL) + { + delete m_vector_view; + } + m_vector_view = new VectorView(m_vector, base_type::total_samples()); + assert(m_vector_view != NULL); + } + + + ~lsf_with_fixed_terms<model_type, true>() override + { + if (m_vector_view != NULL) + { + delete m_vector_view; + } + } + + protected: + Vector m_vector; + VectorView* m_vector_view; + +}; // end class lsf_with_fixed_terms<ModelT, true> + + +} // end namespace detail + + + + +template< typename ModelT, + typename ValueType = typename ModelT::value_type, + bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS > +class least_squeares_fitter +{ +}; + + +template< typename ModelT, typename ValueType > +class least_squeares_fitter<ModelT, ValueType, false> + : public detail::lsf_with_fixed_terms<ModelT> +{ + public: + typedef ModelT model_type; + typedef detail::lsf_with_fixed_terms<model_type> base_type; + typedef typename base_type::parameter_type parameter_type; + typedef typename base_type::value_type value_type; + typedef typename base_type::solution_type solution_type; + + public: + least_squeares_fitter<ModelT, ValueType, false>( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples) + { + } +}; // end class least_squeares_fitter<ModelT, ValueType, true> + + +template< typename ModelT> +class least_squeares_fitter<ModelT, double, true> + : public detail::lsf_with_fixed_terms<ModelT> +{ + public: + typedef ModelT model_type; + typedef detail::lsf_with_fixed_terms<model_type> base_type; + typedef typename base_type::parameter_type parameter_type; + typedef typename base_type::value_type value_type; + typedef typename base_type::solution_type solution_type; + + using base_type::m_vector_view; + //using base_type::result; // VSC legacy support + solution_type& result( std::vector<Point> const& old_sample_values, + std::vector<Point> const& new_sample_values ) + { + return base_type::result(old_sample_values, new_sample_values); + } + + solution_type& result() + { + return base_type::result(); + } + + public: + least_squeares_fitter<ModelT, double, true>( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples) + { + } + + template< typename VectorT > + solution_type& result(VectorT const& sample_values) + { + assert(sample_values.size() == m_vector_view->size()); + Vector sv(sample_values.size()); + for (size_t i = 0; i < sv.size(); ++i) + sv[i] = sample_values[i] - (*m_vector_view)[i]; + return base_type::result(sv); + } + +}; // end class least_squeares_fitter<ModelT, double, true> + + +template< typename ModelT> +class least_squeares_fitter<ModelT, Point, true> + : public detail::lsf_with_fixed_terms<ModelT> +{ + public: + typedef ModelT model_type; + typedef detail::lsf_with_fixed_terms<model_type> base_type; + typedef typename base_type::parameter_type parameter_type; + typedef typename base_type::value_type value_type; + typedef typename base_type::solution_type solution_type; + + using base_type::m_vector_view; + //using base_type::result; // VCS legacy support + solution_type& result( std::vector<Point> const& old_sample_values, + std::vector<Point> const& new_sample_values ) + { + return base_type::result(old_sample_values, new_sample_values); + } + + solution_type& result() + { + return base_type::result(); + } + + + public: + least_squeares_fitter<ModelT, Point, true>( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples) + { + } + + solution_type& result(std::vector<Point> const& sample_values) + { + assert(sample_values.size() == m_vector_view->size()); + NL::Matrix sv(sample_values.size(), 2); + for (size_t i = 0; i < sample_values.size(); ++i) + { + sv(i, X) = sample_values[i][X] - (*m_vector_view)[i]; + sv(i, Y) = sample_values[i][Y] - (*m_vector_view)[i]; + } + return base_type::result(sv); + } + +}; // end class least_squeares_fitter<ModelT, Point, true> + + +} // end namespace NL +} // end namespace Geom + + + +#endif // _NL_FITTING_TOOL_H_ + + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : diff --git a/include/2geom/numeric/linear_system.h b/include/2geom/numeric/linear_system.h new file mode 100644 index 0000000..f793e20 --- /dev/null +++ b/include/2geom/numeric/linear_system.h @@ -0,0 +1,138 @@ +/* + * LinearSystem class wraps some gsl routines for solving linear systems + * + * Authors: + * Marco Cecchetti <mrcekets at gmail.com> + * + * Copyright 2008 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + + +#ifndef _NL_LINEAR_SYSTEM_H_ +#define _NL_LINEAR_SYSTEM_H_ + + +#include <cassert> + +#include <gsl/gsl_linalg.h> + +#include <2geom/numeric/matrix.h> +#include <2geom/numeric/vector.h> + + +namespace Geom { namespace NL { + + +class LinearSystem +{ +public: + LinearSystem(MatrixView & _matrix, VectorView & _vector) + : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns()) + { + } + + LinearSystem(Matrix & _matrix, Vector & _vector) + : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns()) + { + } + + const Vector & LU_solve() + { + assert( matrix().rows() == matrix().columns() + && matrix().rows() == vector().size() ); + int s; + gsl_permutation * p = gsl_permutation_alloc(matrix().rows()); + gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s); + gsl_linalg_LU_solve( matrix().get_gsl_matrix(), + p, + vector().get_gsl_vector(), + m_solution.get_gsl_vector() + ); + gsl_permutation_free(p); + return solution(); + } + + const Vector & SV_solve() + { + assert( matrix().rows() >= matrix().columns() + && matrix().rows() == vector().size() ); + + gsl_matrix* U = matrix().get_gsl_matrix(); + gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns()); + gsl_vector* S = gsl_vector_alloc(matrix().columns()); + gsl_vector* work = gsl_vector_alloc(matrix().columns()); + + gsl_linalg_SV_decomp( U, V, S, work ); + + gsl_vector* b = vector().get_gsl_vector(); + gsl_vector* x = m_solution.get_gsl_vector(); + + gsl_linalg_SV_solve( U, V, S, b, x); + + gsl_matrix_free(V); + gsl_vector_free(S); + gsl_vector_free(work); + + return solution(); + } + + MatrixView & matrix() + { + return m_matrix; + } + + VectorView & vector() + { + return m_vector; + } + + const Vector & solution() const + { + return m_solution; + } + +private: + MatrixView m_matrix; + VectorView m_vector; + Vector m_solution; +}; + + +} } // end namespaces + + +#endif /*_NL_LINEAR_SYSTEM_H_*/ + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : diff --git a/include/2geom/numeric/matrix.h b/include/2geom/numeric/matrix.h new file mode 100644 index 0000000..02851b4 --- /dev/null +++ b/include/2geom/numeric/matrix.h @@ -0,0 +1,603 @@ +/* + * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines; + * "views" mimic the semantic of C++ references: any operation performed + * on a "view" is actually performed on the "viewed object" + * + * Authors: + * Marco Cecchetti <mrcekets at gmail.com> + * + * Copyright 2008 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + + + + +#ifndef _NL_MATRIX_H_ +#define _NL_MATRIX_H_ + +#include <2geom/exception.h> +#include <2geom/numeric/vector.h> + +#include <cassert> +#include <utility> // for std::pair +#include <algorithm> // for std::swap +#include <sstream> +#include <string> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_linalg.h> + + +namespace Geom { namespace NL { + +namespace detail +{ + +class BaseMatrixImpl +{ + public: + virtual ~BaseMatrixImpl() + { + } + + ConstVectorView row_const_view(size_t i) const + { + return ConstVectorView(gsl_matrix_const_row(m_matrix, i)); + } + + ConstVectorView column_const_view(size_t i) const + { + return ConstVectorView(gsl_matrix_const_column(m_matrix, i)); + } + + const double & operator() (size_t i, size_t j) const + { + return *gsl_matrix_const_ptr(m_matrix, i, j); + } + + const gsl_matrix* get_gsl_matrix() const + { + return m_matrix; + } + + bool is_zero() const + { + return gsl_matrix_isnull(m_matrix); + } + + bool is_positive() const + { + for ( unsigned int i = 0; i < rows(); ++i ) + { + for ( unsigned int j = 0; j < columns(); ++j ) + { + if ( (*this)(i,j) <= 0 ) return false; + } + } + return true; + } + + bool is_negative() const + { + for ( unsigned int i = 0; i < rows(); ++i ) + { + for ( unsigned int j = 0; j < columns(); ++j ) + { + if ( (*this)(i,j) >= 0 ) return false; + } + } + return true; + } + + bool is_non_negative() const + { + for ( unsigned int i = 0; i < rows(); ++i ) + { + for ( unsigned int j = 0; j < columns(); ++j ) + { + if ( (*this)(i,j) < 0 ) return false; + } + } + return true; + } + + double max() const + { + return gsl_matrix_max(m_matrix); + } + + double min() const + { + return gsl_matrix_min(m_matrix); + } + + std::pair<size_t, size_t> + max_index() const + { + std::pair<size_t, size_t> indices; + gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second)); + return indices; + } + + std::pair<size_t, size_t> + min_index() const + { + std::pair<size_t, size_t> indices; + gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second)); + return indices; + } + + size_t rows() const + { + return m_rows; + } + + size_t columns() const + { + return m_columns; + } + + std::string str() const; + + protected: + size_t m_rows, m_columns; + gsl_matrix* m_matrix; + +}; // end class BaseMatrixImpl + + +inline +bool operator== (BaseMatrixImpl const& m1, BaseMatrixImpl const& m2) +{ + if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) return false; + + for (size_t i = 0; i < m1.rows(); ++i) + for (size_t j = 0; j < m1.columns(); ++j) + if (m1(i,j) != m2(i,j)) return false; + + return true; +} + +template< class charT > +inline +std::basic_ostream<charT> & +operator<< (std::basic_ostream<charT> & os, const BaseMatrixImpl & _matrix) +{ + if (_matrix.rows() == 0 || _matrix.columns() == 0) return os; + + os << "[[" << _matrix(0,0); + for (size_t j = 1; j < _matrix.columns(); ++j) + { + os << ", " << _matrix(0,j); + } + os << "]"; + + for (size_t i = 1; i < _matrix.rows(); ++i) + { + os << ", [" << _matrix(i,0); + for (size_t j = 1; j < _matrix.columns(); ++j) + { + os << ", " << _matrix(i,j); + } + os << "]"; + } + os << "]"; + return os; +} + +inline +std::string BaseMatrixImpl::str() const +{ + std::ostringstream oss; + oss << (*this); + return oss.str(); +} + + +class MatrixImpl : public BaseMatrixImpl +{ + public: + + typedef BaseMatrixImpl base_type; + + void set_all( double x ) + { + gsl_matrix_set_all(m_matrix, x); + } + + void set_identity() + { + gsl_matrix_set_identity(m_matrix); + } + + using base_type::operator(); // VSC legacy support + const double & operator() (size_t i, size_t j) const + { + return base_type::operator ()(i, j); + } + + double & operator() (size_t i, size_t j) + { + return *gsl_matrix_ptr(m_matrix, i, j); + } + + using base_type::get_gsl_matrix; + + gsl_matrix* get_gsl_matrix() + { + return m_matrix; + } + + VectorView row_view(size_t i) + { + return VectorView(gsl_matrix_row(m_matrix, i)); + } + + VectorView column_view(size_t i) + { + return VectorView(gsl_matrix_column(m_matrix, i)); + } + + void swap_rows(size_t i, size_t j) + { + gsl_matrix_swap_rows(m_matrix, i, j); + } + + void swap_columns(size_t i, size_t j) + { + gsl_matrix_swap_columns(m_matrix, i, j); + } + + MatrixImpl & transpose() + { + assert(columns() == rows()); + gsl_matrix_transpose(m_matrix); + return (*this); + } + + MatrixImpl & scale(double x) + { + gsl_matrix_scale(m_matrix, x); + return (*this); + } + + MatrixImpl & translate(double x) + { + gsl_matrix_add_constant(m_matrix, x); + return (*this); + } + + MatrixImpl & operator+=(base_type const& _matrix) + { + gsl_matrix_add(m_matrix, _matrix.get_gsl_matrix()); + return (*this); + } + + MatrixImpl & operator-=(base_type const& _matrix) + { + gsl_matrix_sub(m_matrix, _matrix.get_gsl_matrix()); + return (*this); + } + +}; // end class MatrixImpl + +} // end namespace detail + + +using detail::operator==; +using detail::operator<<; + + +template <size_t N> +class ConstBaseSymmetricMatrix; + + +class Matrix: public detail::MatrixImpl +{ + public: + typedef detail::MatrixImpl base_type; + + public: + // the matrix is not initialized + Matrix(size_t n1, size_t n2) + { + m_rows = n1; + m_columns = n2; + m_matrix = gsl_matrix_alloc(n1, n2); + } + + Matrix(size_t n1, size_t n2, double x) + { + m_rows = n1; + m_columns = n2; + m_matrix = gsl_matrix_alloc(n1, n2); + gsl_matrix_set_all(m_matrix, x); + } + + Matrix(Matrix const& _matrix) + : base_type() + { + m_rows = _matrix.rows(); + m_columns = _matrix.columns(); + m_matrix = gsl_matrix_alloc(rows(), columns()); + gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix()); + } + + explicit + Matrix(base_type::base_type const& _matrix) + { + m_rows = _matrix.rows(); + m_columns = _matrix.columns(); + m_matrix = gsl_matrix_alloc(rows(), columns()); + gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix()); + } + + template <size_t N> + explicit + Matrix(ConstBaseSymmetricMatrix<N> const& _smatrix) + { + m_rows = N; + m_columns = N; + m_matrix = gsl_matrix_alloc(N, N); + for (size_t i = 0; i < N; ++i) + for (size_t j = 0; j < N ; ++j) + (*gsl_matrix_ptr(m_matrix, i, j)) = _smatrix(i,j); + } + + Matrix & operator=(Matrix const& _matrix) + { + assert( rows() == _matrix.rows() && columns() == _matrix.columns() ); + gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix()); + return *this; + } + + Matrix & operator=(base_type::base_type const& _matrix) + { + assert( rows() == _matrix.rows() && columns() == _matrix.columns() ); + gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix()); + return *this; + } + + template <size_t N> + Matrix & operator=(ConstBaseSymmetricMatrix<N> const& _smatrix) + { + assert (rows() == N && columns() == N); + for (size_t i = 0; i < N; ++i) + for (size_t j = 0; j < N ; ++j) + (*this)(i,j) = _smatrix(i,j); + return *this; + } + + ~Matrix() override + { + gsl_matrix_free(m_matrix); + } + + Matrix & transpose() + { + return static_cast<Matrix &>( base_type::transpose() ); + } + + Matrix & scale(double x) + { + return static_cast<Matrix &>( base_type::scale(x) ); + } + + Matrix & translate(double x) + { + return static_cast<Matrix &>( base_type::translate(x) ); + } + + Matrix & operator+=(base_type::base_type const& _matrix) + { + return static_cast<Matrix &>( base_type::operator+=(_matrix) ); + } + + Matrix & operator-=(base_type::base_type const& _matrix) + { + return static_cast<Matrix &>( base_type::operator-=(_matrix) ); + } + + friend + void swap(Matrix & m1, Matrix & m2); + friend + void swap_any(Matrix & m1, Matrix & m2); + +}; // end class Matrix + + +// warning! this operation invalidates any view of the passed matrix objects +inline +void swap(Matrix & m1, Matrix & m2) +{ + assert(m1.rows() == m2.rows() && m1.columns() == m2.columns()); + using std::swap; + swap(m1.m_matrix, m2.m_matrix); +} + +inline void swap_any(Matrix &m1, Matrix &m2) +{ + using std::swap; + swap(m1.m_matrix, m2.m_matrix); + swap(m1.m_rows, m2.m_rows); + swap(m1.m_columns, m2.m_columns); +} + + + +class ConstMatrixView : public detail::BaseMatrixImpl +{ + public: + typedef detail::BaseMatrixImpl base_type; + + public: + ConstMatrixView(const base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2) + : m_matrix_view( gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2) ) + { + m_rows = n1; + m_columns = n2; + m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) ); + } + + ConstMatrixView(const ConstMatrixView & _matrix) + : base_type(), + m_matrix_view(_matrix.m_matrix_view) + { + m_rows = _matrix.rows(); + m_columns = _matrix.columns(); + m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) ); + } + + ConstMatrixView(const base_type & _matrix) + : m_matrix_view(gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), 0, 0, _matrix.rows(), _matrix.columns())) + { + m_rows = _matrix.rows(); + m_columns = _matrix.columns(); + m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) ); + } + + private: + gsl_matrix_const_view m_matrix_view; + +}; // end class ConstMatrixView + + + + +class MatrixView : public detail::MatrixImpl +{ + public: + typedef detail::MatrixImpl base_type; + + public: + MatrixView(base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2) + { + m_rows = n1; + m_columns = n2; + m_matrix_view + = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2); + m_matrix = &(m_matrix_view.matrix); + } + + MatrixView(const MatrixView & _matrix) + : base_type() + { + m_rows = _matrix.rows(); + m_columns = _matrix.columns(); + m_matrix_view = _matrix.m_matrix_view; + m_matrix = &(m_matrix_view.matrix); + } + + MatrixView(Matrix & _matrix) + { + m_rows = _matrix.rows(); + m_columns = _matrix.columns(); + m_matrix_view + = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), 0, 0, rows(), columns()); + m_matrix = &(m_matrix_view.matrix); + } + + MatrixView & operator=(MatrixView const& _matrix) + { + assert( rows() == _matrix.rows() && columns() == _matrix.columns() ); + gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); + return *this; + } + + MatrixView & operator=(base_type::base_type const& _matrix) + { + assert( rows() == _matrix.rows() && columns() == _matrix.columns() ); + gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix()); + return *this; + } + + MatrixView & transpose() + { + return static_cast<MatrixView &>( base_type::transpose() ); + } + + MatrixView & scale(double x) + { + return static_cast<MatrixView &>( base_type::scale(x) ); + } + + MatrixView & translate(double x) + { + return static_cast<MatrixView &>( base_type::translate(x) ); + } + + MatrixView & operator+=(base_type::base_type const& _matrix) + { + return static_cast<MatrixView &>( base_type::operator+=(_matrix) ); + } + + MatrixView & operator-=(base_type::base_type const& _matrix) + { + return static_cast<MatrixView &>( base_type::operator-=(_matrix) ); + } + + friend + void swap_view(MatrixView & m1, MatrixView & m2); + + private: + gsl_matrix_view m_matrix_view; + +}; // end class MatrixView + + +inline +void swap_view(MatrixView & m1, MatrixView & m2) +{ + assert(m1.rows() == m2.rows() && m1.columns() == m2.columns()); + using std::swap; + swap(m1.m_matrix_view, m2.m_matrix_view); +} + +Vector operator*( detail::BaseMatrixImpl const& A, + detail::BaseVectorImpl const& v ); + +Matrix operator*( detail::BaseMatrixImpl const& A, + detail::BaseMatrixImpl const& B ); + +Matrix pseudo_inverse(detail::BaseMatrixImpl const& A); + +double trace (detail::BaseMatrixImpl const& A); + +double det (detail::BaseMatrixImpl const& A); + +} } // end namespaces + +#endif /*_NL_MATRIX_H_*/ + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : diff --git a/include/2geom/numeric/symmetric-matrix-fs-operation.h b/include/2geom/numeric/symmetric-matrix-fs-operation.h new file mode 100644 index 0000000..c5aaa72 --- /dev/null +++ b/include/2geom/numeric/symmetric-matrix-fs-operation.h @@ -0,0 +1,102 @@ +/* + * SymmetricMatrix basic operation + * + * Authors: + * Marco Cecchetti <mrcekets at gmail.com> + * + * Copyright 2009 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + +#ifndef _NL_SYMMETRIC_MATRIX_FS_OPERATION_H_ +#define _NL_SYMMETRIC_MATRIX_FS_OPERATION_H_ + + +#include <2geom/numeric/symmetric-matrix-fs.h> +#include <2geom/numeric/symmetric-matrix-fs-trace.h> + + + + +namespace Geom { namespace NL { + +template <size_t N> +SymmetricMatrix<N> adj(const ConstBaseSymmetricMatrix<N> & S); + +template <> +inline +SymmetricMatrix<2> adj(const ConstBaseSymmetricMatrix<2> & S) +{ + SymmetricMatrix<2> result; + result.get<0,0>() = S.get<1,1>(); + result.get<1,0>() = -S.get<1,0>(); + result.get<1,1>() = S.get<0,0>(); + return result; +} + +template <> +inline +SymmetricMatrix<3> adj(const ConstBaseSymmetricMatrix<3> & S) +{ + SymmetricMatrix<3> result; + + result.get<0,0>() = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>(); + result.get<1,0>() = S.get<0,2>() * S.get<2,1>() - S.get<0,1>() * S.get<2,2>(); + result.get<1,1>() = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>(); + result.get<2,0>() = S.get<0,1>() * S.get<1,2>() - S.get<0,2>() * S.get<1,1>(); + result.get<2,1>() = S.get<0,2>() * S.get<1,0>() - S.get<0,0>() * S.get<1,2>(); + result.get<2,2>() = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>(); + return result; +} + +template <size_t N> +inline +SymmetricMatrix<N> inverse(const ConstBaseSymmetricMatrix<N> & S) +{ + SymmetricMatrix<N> result = adj(S); + double d = det(S); + assert (d != 0); + result.scale (1/d); + return result; +} + +} /* end namespace NL*/ } /* end namespace Geom*/ + + +#endif // _NL_SYMMETRIC_MATRIX_FS_OPERATION_H_ + + + + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : diff --git a/include/2geom/numeric/symmetric-matrix-fs-trace.h b/include/2geom/numeric/symmetric-matrix-fs-trace.h new file mode 100644 index 0000000..0e7a28c --- /dev/null +++ b/include/2geom/numeric/symmetric-matrix-fs-trace.h @@ -0,0 +1,427 @@ +/* + * SymmetricMatrix trace + * + * Authors: + * Marco Cecchetti <mrcekets at gmail.com> + * + * Copyright 2009 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + + +#ifndef _NL_TRACE_H_ +#define _NL_TRACE_H_ + + +#include <2geom/numeric/matrix.h> +#include <2geom/numeric/symmetric-matrix-fs.h> + + + + + +namespace Geom { namespace NL { + + +namespace detail +{ + +/* + * helper routines + */ + +inline +int sgn_prod (int x, int y) +{ + if (x == 0 || y == 0) return 0; + if (x == y) return 1; + return -1; +} + +inline +bool abs_less (double x, double y) +{ + return (std::fabs(x) < std::fabs(y)); +} + + +/* + * trace K-th of symmetric matrix S of order N + */ +template <size_t K, size_t N> +struct trace +{ + static double evaluate(const ConstBaseSymmetricMatrix<N> &S); +}; + +template <size_t N> +struct trace<1,N> +{ + static + double evaluate (const ConstBaseSymmetricMatrix<N> & S) + { + double t = 0; + for (size_t i = 0; i < N; ++i) + { + t += S(i,i); + } + return t; + } +}; + +template <size_t N> +struct trace<N,N> +{ + static + double evaluate (const ConstBaseSymmetricMatrix<N> & S) + { + Matrix M(S); + return det(M); + } +}; + +/* + * trace for symmetric matrix of order 2 + */ +template <> +struct trace<1,2> +{ + static + double evaluate (const ConstBaseSymmetricMatrix<2> & S) + { + return (S.get<0,0>() + S.get<1,1>()); + } +}; + +template <> +struct trace<2,2> +{ + static + double evaluate (const ConstBaseSymmetricMatrix<2> & S) + { + return (S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>()); + } +}; + + +/* + * trace for symmetric matrix of order 3 + */ +template <> +struct trace<1,3> +{ + static + double evaluate (const ConstBaseSymmetricMatrix<3> & S) + { + return (S.get<0,0>() + S.get<1,1>() + S.get<2,2>()); + } +}; + +template <> +struct trace<2,3> +{ + static + double evaluate (const ConstBaseSymmetricMatrix<3> & S) + { + double a00 = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>(); + double a11 = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>(); + double a22 = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>(); + return (a00 + a11 + a22); + } +}; + +template <> +struct trace<3,3> +{ + static + double evaluate (const ConstBaseSymmetricMatrix<3> & S) + { + double d = S.get<0,0>() * S.get<1,1>() * S.get<2,2>(); + d += (2 * S.get<1,0>() * S.get<2,0>() * S.get<2,1>()); + d -= (S.get<0,0>() * S.get<2,1>() * S.get<2,1>()); + d -= (S.get<1,1>() * S.get<2,0>() * S.get<2,0>()); + d -= (S.get<2,2>() * S.get<1,0>() * S.get<1,0>()); + return d; + } +}; + + +/* + * sign of trace K-th + */ +template <size_t K, size_t N> +struct trace_sgn +{ + static + int evaluate (const ConstBaseSymmetricMatrix<N> & S) + { + double d = trace<K, N>::evaluate(S); + return sgn(d); + } +}; + + +/* + * sign of trace for symmetric matrix of order 2 + */ +template <> +struct trace_sgn<2,2> +{ + static + int evaluate (const ConstBaseSymmetricMatrix<2> & S) + { + double m00 = S.get<0,0>(); + double m10 = S.get<1,0>(); + double m11 = S.get<1,1>(); + + int sm00 = sgn (m00); + int sm10 = sgn (m10); + int sm11 = sgn (m11); + + if (sm10 == 0) + { + return sgn_prod (sm00, sm11); + } + else + { + int sm00m11 = sgn_prod (sm00, sm11); + if (sm00m11 == 1) + { + int e00, e10, e11; + double f00 = std::frexp (m00, &e00); + double f10 = std::frexp (m10, &e10); + double f11 = std::frexp (m11, &e11); + + int e0011 = e00 + e11; + int e1010 = e10 << 1; + int ed = e0011 - e1010; + + if (ed > 1) + { + return 1; + } + else if (ed < -1) + { + return -1; + } + else + { + double d = std::ldexp (f00 * f11, ed) - f10 * f10; + //std::cout << "trace_sgn<2,2>: det = " << d << std::endl; + double eps = std::ldexp (1, -50); + if (std::fabs(d) < eps) return 0; + return sgn (d); + } + } + return -1; + } + } +}; + + +/* + * sign of trace for symmetric matrix of order 3 + */ +template <> +struct trace_sgn<2,3> +{ + static + int evaluate (const ConstBaseSymmetricMatrix<3> & S) + { + double eps = std::ldexp (1, -50); + double t[6]; + + t[0] = S.get<1,1>() * S.get<2,2>(); + t[1] = - S.get<1,2>() * S.get<2,1>(); + t[2] = S.get<0,0>() * S.get<2,2>(); + t[3] = - S.get<0,2>() * S.get<2,0>(); + t[4] = S.get<0,0>() * S.get<1,1>(); + t[5] = - S.get<0,1>() * S.get<1,0>(); + + + double* maxp = std::max_element (t, t+6, abs_less); + int em; + std::frexp(*maxp, &em); + double d = 0; + for (double i : t) + { + d += i; + } + double r = std::fabs (std::ldexp (d, -em)); // relative error + //std::cout << "trace_sgn<2,3>: d = " << d << std::endl; + //std::cout << "trace_sgn<2,3>: r = " << r << std::endl; + if (r < eps) return 0; + if (d > 0) return 1; + return -1; + } +}; + +template <> +struct trace_sgn<3,3> +{ + static + int evaluate (const ConstBaseSymmetricMatrix<3> & S) + { + + double eps = std::ldexp (1, -48); + double t[5]; + + t[0] = S.get<0,0>() * S.get<1,1>() * S.get<2,2>(); + t[1] = 2 * S.get<1,0>() * S.get<2,0>() * S.get<2,1>(); + t[2] = -(S.get<0,0>() * S.get<2,1>() * S.get<2,1>()); + t[3] = -(S.get<1,1>() * S.get<2,0>() * S.get<2,0>()); + t[4] = -(S.get<2,2>() * S.get<1,0>() * S.get<1,0>()); + + double* maxp = std::max_element (t, t+5, abs_less); + int em; + std::frexp(*maxp, &em); + double d = 0; + for (double i : t) + { + d += i; + } + //std::cout << "trace_sgn<3,3>: d = " << d << std::endl; + double r = std::fabs (std::ldexp (d, -em)); // relative error + //std::cout << "trace_sgn<3,3>: r = " << r << std::endl; + + if (r < eps) return 0; + if (d > 0) return 1; + return -1; + } +}; // end struct trace_sgn<3,3> + +} // end namespace detail + + +template <size_t K, size_t N> +inline +double trace (const ConstBaseSymmetricMatrix<N> & _matrix) +{ + return detail::trace<K, N>::evaluate(_matrix); +} + +template <size_t N> +inline +double trace (const ConstBaseSymmetricMatrix<N> & _matrix) +{ + return detail::trace<1, N>::evaluate(_matrix); +} + +template <size_t N> +inline +double det (const ConstBaseSymmetricMatrix<N> & _matrix) +{ + return detail::trace<N, N>::evaluate(_matrix); +} + + +template <size_t K, size_t N> +inline +int trace_sgn (const ConstBaseSymmetricMatrix<N> & _matrix) +{ + return detail::trace_sgn<K, N>::evaluate(_matrix); +} + +template <size_t N> +inline +int trace_sgn (const ConstBaseSymmetricMatrix<N> & _matrix) +{ + return detail::trace_sgn<1, N>::evaluate(_matrix); +} + +template <size_t N> +inline +int det_sgn (const ConstBaseSymmetricMatrix<N> & _matrix) +{ + return detail::trace_sgn<N, N>::evaluate(_matrix); +} + +/* +template <size_t N> +inline +size_t rank (const ConstBaseSymmetricMatrix<N> & S) +{ + THROW_NOTIMPLEMENTED(); + return 0; +} + +template <> +inline +size_t rank<2> (const ConstBaseSymmetricMatrix<2> & S) +{ + if (S.is_zero()) return 0; + double d = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>(); + if (d != 0) return 2; + return 1; +} + +template <> +inline +size_t rank<3> (const ConstBaseSymmetricMatrix<3> & S) +{ + if (S.is_zero()) return 0; + + double a20 = S.get<0,1>() * S.get<1,2>() - S.get<0,2>() * S.get<1,1>(); + double a21 = S.get<0,2>() * S.get<1,0>() - S.get<0,0>() * S.get<1,2>(); + double a22 = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>(); + double d = a20 * S.get<2,0>() + a21 * S.get<2,1>() + a22 * S.get<2,2>(); + + if (d != 0) return 3; + + if (a20 != 0 || a21 != 0 || a22 != 0) return 2; + + double a00 = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>(); + if (a00 != 0) return 2; + + double a10 = S.get<0,2>() * S.get<2,1>() - S.get<0,1>() * S.get<2,2>(); + if (a10 != 0) return 2; + + double a11 = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>(); + if (a11 != 0) return 2; + + return 1; +} +*/ + +} /* end namespace NL*/ } /* end namespace Geom*/ + + + + +#endif // _NL_TRACE_H_ + + + + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : diff --git a/include/2geom/numeric/symmetric-matrix-fs.h b/include/2geom/numeric/symmetric-matrix-fs.h new file mode 100644 index 0000000..2fadd69 --- /dev/null +++ b/include/2geom/numeric/symmetric-matrix-fs.h @@ -0,0 +1,733 @@ +/* + * SymmetricMatrix, ConstSymmetricMatrixView, SymmetricMatrixView template + * classes implement fixed size symmetric matrix; "views" mimic the semantic + * of C++ references: any operation performed on a "view" is actually performed + * on the "viewed object" + * + * Authors: + * Marco Cecchetti <mrcekets at gmail.com> + * + * Copyright 2009 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + + +#ifndef _NL_SYMMETRIC_MATRIX_FS_H_ +#define _NL_SYMMETRIC_MATRIX_FS_H_ + + +#include <2geom/numeric/vector.h> +#include <2geom/numeric/matrix.h> +#include <2geom/utils.h> +#include <2geom/exception.h> + +#include <boost/static_assert.hpp> + +#include <cassert> +#include <utility> // for std::pair +#include <algorithm> // for std::swap, std::copy +#include <sstream> +#include <string> + + + +namespace Geom { namespace NL { + + +namespace detail +{ + +template <size_t I, size_t J, bool B = (I < J)> +struct index +{ + static const size_t K = index<J, I>::K; +}; + +template <size_t I, size_t J> +struct index<I, J, false> +{ + static const size_t K = (((I+1) * I) >> 1) + J; +}; + +} // end namespace detail + + + + +template <size_t N> +class ConstBaseSymmetricMatrix; + +template <size_t N> +class BaseSymmetricMatrix; + +template <size_t N> +class SymmetricMatrix; + +template <size_t N> +class ConstSymmetricMatrixView; + +template <size_t N> +class SymmetricMatrixView; + + + +// declaration needed for friend clause +template <size_t N> +bool operator== (ConstBaseSymmetricMatrix<N> const& _smatrix1, + ConstBaseSymmetricMatrix<N> const& _smatrix2); + + + + +template <size_t N> +class ConstBaseSymmetricMatrix +{ + public: + const static size_t DIM = N; + const static size_t DATA_SIZE = ((DIM+1) * DIM) / 2; + + public: + + ConstBaseSymmetricMatrix (VectorView const& _data) + : m_data(_data) + { + } + + double operator() (size_t i, size_t j) const + { + return m_data[get_index(i,j)]; + } + + template <size_t I, size_t J> + double get() const + { + BOOST_STATIC_ASSERT ((I < N && J < N)); + return m_data[detail::index<I, J>::K]; + } + + + size_t rows() const + { + return DIM; + } + + size_t columns() const + { + return DIM; + } + + bool is_zero() const + { + return m_data.is_zero(); + } + + bool is_positive() const + { + return m_data.is_positive(); + } + + bool is_negative() const + { + return m_data.is_negative(); + } + + bool is_non_negative() const + { + return m_data.is_non_negative(); + } + + double min() const + { + return m_data.min(); + } + + double max() const + { + return m_data.max(); + } + + std::pair<size_t, size_t> + min_index() const + { + std::pair<size_t, size_t> indices(0,0); + double min_value = m_data[0]; + for (size_t i = 1; i < DIM; ++i) + { + for (size_t j = 0; j <= i; ++j) + { + if (min_value > (*this)(i,j)) + { + min_value = (*this)(i,j); + indices.first = i; + indices.second = j; + } + } + } + return indices; + } + + std::pair<size_t, size_t> + max_index() const + { + std::pair<size_t, size_t> indices(0,0); + double max_value = m_data[0]; + for (size_t i = 1; i < DIM; ++i) + { + for (size_t j = 0; j <= i; ++j) + { + if (max_value < (*this)(i,j)) + { + max_value = (*this)(i,j); + indices.first = i; + indices.second = j; + } + } + } + return indices; + } + + size_t min_on_row_index (size_t i) const + { + size_t idx = 0; + double min_value = (*this)(i,0); + for (size_t j = 1; j < DIM; ++j) + { + if (min_value > (*this)(i,j)) + { + min_value = (*this)(i,j); + idx = j; + } + } + return idx; + } + + size_t max_on_row_index (size_t i) const + { + size_t idx = 0; + double max_value = (*this)(i,0); + for (size_t j = 1; j < DIM; ++j) + { + if (max_value < (*this)(i,j)) + { + max_value = (*this)(i,j); + idx = j; + } + } + return idx; + } + + size_t min_on_column_index (size_t j) const + { + return min_on_row_index(j); + } + + size_t max_on_column_index (size_t j) const + { + return max_on_row_index(j); + } + + size_t min_on_diag_index () const + { + size_t idx = 0; + double min_value = (*this)(0,0); + for (size_t i = 1; i < DIM; ++i) + { + if (min_value > (*this)(i,i)) + { + min_value = (*this)(i,i); + idx = i; + } + } + return idx; + } + + size_t max_on_diag_index () const + { + size_t idx = 0; + double max_value = (*this)(0,0); + for (size_t i = 1; i < DIM; ++i) + { + if (max_value < (*this)(i,i)) + { + max_value = (*this)(i,i); + idx = i; + } + } + return idx; + } + + std::string str() const; + + ConstSymmetricMatrixView<N-1> main_minor_const_view() const; + + SymmetricMatrix<N> operator- () const; + + Vector operator* (ConstVectorView _vector) const + { + assert (_vector.size() == DIM); + Vector result(DIM, 0.0); + + for (size_t i = 0; i < DIM; ++i) + { + for (size_t j = 0; j < DIM; ++j) + { + result[i] += (*this)(i,j) * _vector[j]; + } + } + return result; + } + + protected: + static size_t get_index (size_t i, size_t j) + { + if (i < j) return get_index (j, i); + size_t k = (i+1) * i; + k >>= 1; + k += j; + return k; + } + + protected: + ConstVectorView get_data() const + { + return m_data; + } + + friend + bool operator==<N> (ConstBaseSymmetricMatrix const& _smatrix1, + ConstBaseSymmetricMatrix const& _smatrix2); + + protected: + VectorView m_data; + +}; //end ConstBaseSymmetricMatrix + + +template <size_t N> +class BaseSymmetricMatrix : public ConstBaseSymmetricMatrix<N> +{ + public: + typedef ConstBaseSymmetricMatrix<N> base_type; + + + public: + + BaseSymmetricMatrix (VectorView const& _data) + : base_type(_data) + { + } + + double operator() (size_t i, size_t j) const + { + return m_data[base_type::get_index(i,j)]; + } + + double& operator() (size_t i, size_t j) + { + return m_data[base_type::get_index(i,j)]; + } + + template <size_t I, size_t J> + double& get() + { + BOOST_STATIC_ASSERT ((I < N && J < N)); + return m_data[detail::index<I, J>::K]; + } + + void set_all (double x) + { + m_data.set_all(x); + } + + SymmetricMatrixView<N-1> main_minor_view(); + + BaseSymmetricMatrix& transpose() const + { + return (*this); + } + + BaseSymmetricMatrix& translate (double c) + { + m_data.translate(c); + return (*this); + } + + BaseSymmetricMatrix& scale (double c) + { + m_data.scale(c); + return (*this); + } + + BaseSymmetricMatrix& operator+= (base_type const& _smatrix) + { + m_data += (static_cast<const BaseSymmetricMatrix &>(_smatrix).m_data); + return (*this); + } + + BaseSymmetricMatrix& operator-= (base_type const& _smatrix) + { + m_data -= (static_cast<const BaseSymmetricMatrix &>(_smatrix).m_data); + return (*this); + } + + using base_type::DIM; + using base_type::DATA_SIZE; + using base_type::m_data; + using base_type::operator-; + using base_type::operator*; + +}; //end BaseSymmetricMatrix + + +template <size_t N> +class SymmetricMatrix : public BaseSymmetricMatrix<N> +{ + public: + typedef BaseSymmetricMatrix<N> base_type; + typedef typename base_type::base_type base_base_type; + + using base_type::DIM; + using base_type::DATA_SIZE; + using base_type::m_data; + + public: + SymmetricMatrix () + : base_type (VectorView(m_adata, DATA_SIZE)) + { + } + + explicit + SymmetricMatrix (ConstVectorView _data) + : base_type (VectorView(m_adata, DATA_SIZE)) + { + assert (_data.size() == DATA_SIZE); + m_data = _data; + } + + explicit + SymmetricMatrix (const double _data[DATA_SIZE]) + : base_type (VectorView(m_adata, DATA_SIZE)) + { + std::copy (_data, _data + DATA_SIZE, m_adata); + } + + SymmetricMatrix (SymmetricMatrix const& _smatrix) + : base_type (VectorView(m_adata, DATA_SIZE)) + { + m_data = _smatrix.m_data; + } + + explicit + SymmetricMatrix (base_base_type const& _smatrix) + : base_type (VectorView(m_adata, DATA_SIZE)) + { + m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data; + } + + explicit + SymmetricMatrix (ConstMatrixView const& _matrix) + : base_type (VectorView(m_adata, DATA_SIZE)) + { + assert (_matrix.rows() == N && _matrix.columns() == N); + for (size_t i = 0; i < N; ++i) + for (size_t j = 0; j <= i ; ++j) + (*this)(i,j) = _matrix(i,j); + } + + SymmetricMatrix& operator= (SymmetricMatrix const& _smatrix) + { + m_data = _smatrix.m_data; + return (*this); + } + + SymmetricMatrix& operator= (base_base_type const& _smatrix) + { + + m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data; + return (*this); + } + + SymmetricMatrix& operator= (ConstMatrixView const& _matrix) + { + assert (_matrix.rows() == N && _matrix.columns() == N); + for (size_t i = 0; i < N; ++i) + for (size_t j = 0; j <= i ; ++j) + (*this)(i,j) = _matrix(i,j); + + return (*this); + } + + // needed for accessing m_adata + friend class ConstSymmetricMatrixView<DIM>; + friend class SymmetricMatrixView<DIM>; + private: + double m_adata[DATA_SIZE]; +}; //end SymmetricMatrix + + +template <size_t N> +class ConstSymmetricMatrixView : public ConstBaseSymmetricMatrix<N> +{ + public: + typedef ConstBaseSymmetricMatrix<N> base_type; + + using base_type::DIM; + using base_type::DATA_SIZE; + using base_type::m_data; + + + public: + + explicit + ConstSymmetricMatrixView (ConstVectorView _data) + : base_type (const_vector_view_cast(_data)) + { + assert (_data.size() == DATA_SIZE); + } + + explicit + ConstSymmetricMatrixView (const double _data[DATA_SIZE]) + : base_type (const_vector_view_cast (ConstVectorView (_data, DATA_SIZE))) + { + } + + ConstSymmetricMatrixView (const ConstSymmetricMatrixView & _smatrix) + : base_type (_smatrix.m_data) + { + } + + ConstSymmetricMatrixView (const base_type & _smatrix) + : base_type (static_cast<const ConstSymmetricMatrixView &>(_smatrix).m_data) + { + } + +}; //end SymmetricMatrix + + +// declaration needed for friend clause +template <size_t N> +void swap_view(SymmetricMatrixView<N> & m1, SymmetricMatrixView<N> & m2); + + +template <size_t N> +class SymmetricMatrixView : public BaseSymmetricMatrix<N> +{ + public: + typedef BaseSymmetricMatrix<N> base_type; + typedef typename base_type::base_type base_base_type; + + using base_type::DIM; + using base_type::DATA_SIZE; + using base_type::m_data; + + public: + + explicit + SymmetricMatrixView (VectorView _data) + : base_type (_data) + { + assert (_data.size() == DATA_SIZE); + } + + explicit + SymmetricMatrixView (double _data[DATA_SIZE]) + : base_type (VectorView (_data, DATA_SIZE)) + { + } + + SymmetricMatrixView (const SymmetricMatrixView & _smatrix) + : base_type (_smatrix.m_data) + { + } + + SymmetricMatrixView (SymmetricMatrix<DIM> & _smatrix) + : base_type (VectorView (_smatrix.m_adata, DATA_SIZE)) + { + } + + SymmetricMatrixView& operator= (const SymmetricMatrixView & _smatrix) + { + m_data = _smatrix.m_data; + return (*this); + } + + SymmetricMatrixView& operator= (const base_base_type & _smatrix) + { + m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data; + return (*this); + } + + friend + void swap_view<N>(SymmetricMatrixView & m1, SymmetricMatrixView & m2); + +}; //end SymmetricMatrix + + + + +/* + * class ConstBaseSymmetricMatrix methods + */ + +template <size_t N> +inline +std::string ConstBaseSymmetricMatrix<N>::str() const +{ + std::ostringstream oss; + oss << (*this); + return oss.str(); +} + +template <size_t N> +inline +ConstSymmetricMatrixView<N-1> +ConstBaseSymmetricMatrix<N>::main_minor_const_view() const +{ + ConstVectorView data(m_data.get_gsl_vector()->data, DATA_SIZE - DIM); + ConstSymmetricMatrixView<N-1> mm(data); + return mm; +} + +template <size_t N> +inline +SymmetricMatrix<N> ConstBaseSymmetricMatrix<N>::operator- () const +{ + SymmetricMatrix<N> result; + for (size_t i = 0; i < DATA_SIZE; ++i) + { + result.m_data[i] = -m_data[i]; + } + return result; +} + + +/* + * class ConstBaseSymmetricMatrix friend free functions + */ + +template <size_t N> +inline +bool operator== (ConstBaseSymmetricMatrix<N> const& _smatrix1, + ConstBaseSymmetricMatrix<N> const& _smatrix2) +{ + return (_smatrix1.m_data == _smatrix2.m_data); +} + +/* + * class ConstBaseSymmetricMatrix related free functions + */ + +template< size_t N, class charT > +inline +std::basic_ostream<charT> & +operator<< (std::basic_ostream<charT> & os, + const ConstBaseSymmetricMatrix<N> & _matrix) +{ + os << "[[" << _matrix(0,0); + for (size_t j = 1; j < N; ++j) + { + os << ", " << _matrix(0,j); + } + os << "]"; + for (size_t i = 1; i < N; ++i) + { + os << "\n [" << _matrix(i,0); + for (size_t j = 1; j < N; ++j) + { + os << ", " << _matrix(i,j); + } + os << "]"; + } + os << "]"; + return os; +} + + +/* + * class ConstBaseSymmetricMatrix specialized methods + */ + +template<> +inline +size_t ConstBaseSymmetricMatrix<2>::get_index (size_t i, size_t j) +{ + return (i+j); +} + +template<> +inline +size_t ConstBaseSymmetricMatrix<3>::get_index (size_t i, size_t j) +{ + size_t k = i + j; + if (i == 2 || j == 2) ++k; + return k; +} + + +/* + * class BaseSymmetricMatrix methods + */ + +template <size_t N> +inline +SymmetricMatrixView<N-1> BaseSymmetricMatrix<N>::main_minor_view() +{ + VectorView data(m_data.get_gsl_vector()->data, DATA_SIZE - DIM); + SymmetricMatrixView<N-1> mm(data); + return mm; +} + + +/* + * class SymmetricMatrixView friend free functions + */ + +template <size_t N> +inline +void swap_view(SymmetricMatrixView<N> & m1, SymmetricMatrixView<N> & m2) +{ + swap_view(m1.m_data, m2.m_data); +} + +} /* end namespace NL*/ } /* end namespace Geom*/ + + + + +#endif // _NL_SYMMETRIC_MATRIX_FS_H_ + + + + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : diff --git a/include/2geom/numeric/vector.h b/include/2geom/numeric/vector.h new file mode 100644 index 0000000..b66115b --- /dev/null +++ b/include/2geom/numeric/vector.h @@ -0,0 +1,594 @@ +/* + * Vector, VectorView, ConstVectorView classes wrap the gsl vector routines; + * "views" mimic the semantic of C++ references: any operation performed + * on a "view" is actually performed on the "viewed object" + * + * Authors: + * Marco Cecchetti <mrcekets at gmail.com> + * + * Copyright 2008 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + + + + +#ifndef _NL_VECTOR_H_ +#define _NL_VECTOR_H_ + +#include <cassert> +#include <algorithm> // for std::swap +#include <vector> +#include <sstream> +#include <string> + + +#include <gsl/gsl_vector.h> +#include <gsl/gsl_blas.h> + + +namespace Geom { namespace NL { + +namespace detail +{ + +class BaseVectorImpl +{ + public: + double const& operator[](size_t i) const + { + return *gsl_vector_const_ptr(m_vector, i); + } + + const gsl_vector* get_gsl_vector() const + { + return m_vector; + } + bool is_zero() const + { + return gsl_vector_isnull(m_vector); + } + + bool is_positive() const + { + for ( size_t i = 0; i < size(); ++i ) + { + if ( (*this)[i] <= 0 ) return false; + } + return true; + } + + bool is_negative() const + { + for ( size_t i = 0; i < size(); ++i ) + { + if ( (*this)[i] >= 0 ) return false; + } + return true; + } + + bool is_non_negative() const + { + for ( size_t i = 0; i < size(); ++i ) + { + if ( (*this)[i] < 0 ) return false; + } + return true; + } + + double max() const + { + return gsl_vector_max(m_vector); + } + + double min() const + { + return gsl_vector_min(m_vector); + } + + size_t max_index() const + { + return gsl_vector_max_index(m_vector); + } + + size_t min_index() const + { + return gsl_vector_min_index(m_vector); + } + + size_t size() const + { + return m_size; + } + + std::string str() const; + + virtual ~BaseVectorImpl() + { + } + + protected: + size_t m_size; + gsl_vector* m_vector; + +}; // end class BaseVectorImpl + + +inline +bool operator== (BaseVectorImpl const& v1, BaseVectorImpl const& v2) +{ + if (v1.size() != v2.size()) return false; + + for (size_t i = 0; i < v1.size(); ++i) + { + if (v1[i] != v2[i]) return false; + } + return true; +} + +template< class charT > +inline +std::basic_ostream<charT> & +operator<< (std::basic_ostream<charT> & os, const BaseVectorImpl & _vector) +{ + if (_vector.size() == 0 ) return os; + os << "[" << _vector[0]; + for (unsigned int i = 1; i < _vector.size(); ++i) + { + os << ", " << _vector[i]; + } + os << "]"; + return os; +} + +inline +std::string BaseVectorImpl::str() const +{ + std::ostringstream oss; + oss << (*this); + return oss.str(); +} + +inline +double dot(BaseVectorImpl const& v1, BaseVectorImpl const& v2) +{ + double result; + gsl_blas_ddot(v1.get_gsl_vector(), v2.get_gsl_vector(), &result); + return result; +} + + +class VectorImpl : public BaseVectorImpl +{ + public: + typedef BaseVectorImpl base_type; + + public: + void set_all(double x) + { + gsl_vector_set_all(m_vector, x); + } + + void set_basis(size_t i) + { + gsl_vector_set_basis(m_vector, i); + } + + using base_type::operator[]; + + double & operator[](size_t i) + { + return *gsl_vector_ptr(m_vector, i); + } + + using base_type::get_gsl_vector; + + gsl_vector* get_gsl_vector() + { + return m_vector; + } + + void swap_elements(size_t i, size_t j) + { + gsl_vector_swap_elements(m_vector, i, j); + } + + void reverse() + { + gsl_vector_reverse(m_vector); + } + + VectorImpl & scale(double x) + { + gsl_vector_scale(m_vector, x); + return (*this); + } + + VectorImpl & translate(double x) + { + gsl_vector_add_constant(m_vector, x); + return (*this); + } + + VectorImpl & operator+=(base_type const& _vector) + { + gsl_vector_add(m_vector, _vector.get_gsl_vector()); + return (*this); + } + + VectorImpl & operator-=(base_type const& _vector) + { + gsl_vector_sub(m_vector, _vector.get_gsl_vector()); + return (*this); + } + +}; // end class VectorImpl + +} // end namespace detail + + +using detail::operator==; +using detail::operator<<; + +class Vector : public detail::VectorImpl +{ + public: + typedef detail::VectorImpl base_type; + + public: + Vector(size_t n) + { + m_size = n; + m_vector = gsl_vector_alloc(n); + } + + Vector(size_t n, double x) + { + m_size = n; + m_vector = gsl_vector_alloc(n); + gsl_vector_set_all(m_vector, x); + } + + // create a vector with n elements all set to zero + // but the i-th that is set to 1 + Vector(size_t n, size_t i) + { + m_size = n; + m_vector = gsl_vector_alloc(n); + gsl_vector_set_basis(m_vector, i); + } + + Vector(Vector const& _vector) + : base_type() + { + m_size = _vector.size(); + m_vector = gsl_vector_alloc(size()); + gsl_vector_memcpy(m_vector, _vector.m_vector); + } + + explicit + Vector(base_type::base_type const& _vector) + { + m_size = _vector.size(); + m_vector = gsl_vector_alloc(size()); + gsl_vector_memcpy(m_vector, _vector.get_gsl_vector()); + } + + ~Vector() override + { + gsl_vector_free(m_vector); + } + + + Vector & operator=(Vector const& _vector) + { + assert( size() == _vector.size() ); + gsl_vector_memcpy(m_vector, _vector.m_vector); + return (*this); + } + + Vector & operator=(base_type::base_type const& _vector) + { + assert( size() == _vector.size() ); + gsl_vector_memcpy(m_vector, _vector.get_gsl_vector()); + return (*this); + } + + Vector & scale(double x) + { + return static_cast<Vector&>( base_type::scale(x) ); + } + + Vector & translate(double x) + { + return static_cast<Vector&>( base_type::translate(x) ); + } + + Vector & operator+=(base_type::base_type const& _vector) + { + return static_cast<Vector&>( base_type::operator+=(_vector) ); + } + + Vector & operator-=(base_type::base_type const& _vector) + { + return static_cast<Vector&>( base_type::operator-=(_vector) ); + } + + friend + void swap(Vector & v1, Vector & v2); + friend + void swap_any(Vector & v1, Vector & v2); + +}; // end class Vector + + +// warning! these operations invalidate any view of the passed vector objects +inline +void swap(Vector & v1, Vector & v2) +{ + assert(v1.size() == v2.size()); + using std::swap; + swap(v1.m_vector, v2.m_vector); +} + +inline +void swap_any(Vector & v1, Vector & v2) +{ + using std::swap; + swap(v1.m_vector, v2.m_vector); + swap(v1.m_size, v2.m_size); +} + + +class ConstVectorView : public detail::BaseVectorImpl +{ + public: + typedef detail::BaseVectorImpl base_type; + + public: + ConstVectorView(const base_type & _vector, size_t n, size_t offset = 0) + : m_vector_view( gsl_vector_const_subvector(_vector.get_gsl_vector(), offset, n) ) + { + m_size = n; + m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) ); + } + + ConstVectorView(const base_type & _vector, size_t n, size_t offset , size_t stride) + : m_vector_view( gsl_vector_const_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n) ) + { + m_size = n; + m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) ); + } + + ConstVectorView(const double* _vector, size_t n, size_t offset = 0) + : m_vector_view( gsl_vector_const_view_array(_vector + offset, n) ) + { + m_size = n; + m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) ); + } + + ConstVectorView(const double* _vector, size_t n, size_t offset, size_t stride) + : m_vector_view( gsl_vector_const_view_array_with_stride(_vector + offset, stride, n) ) + { + m_size = n; + m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) ); + } + + explicit + ConstVectorView(gsl_vector_const_view _gsl_vector_view) + : m_vector_view(_gsl_vector_view) + { + m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) ); + m_size = m_vector->size; + } + + explicit + ConstVectorView(const std::vector<double>& _vector) + : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) ) + { + m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) ); + m_size = _vector.size(); + } + + ConstVectorView(const ConstVectorView & _vector) + : base_type(), + m_vector_view(_vector.m_vector_view) + { + m_size = _vector.size(); + m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) ); + } + + ConstVectorView(const base_type & _vector) + : m_vector_view(gsl_vector_const_subvector(_vector.get_gsl_vector(), 0, _vector.size())) + { + m_size = _vector.size(); + m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) ); + } + + private: + gsl_vector_const_view m_vector_view; + +}; // end class ConstVectorView + + + + +class VectorView : public detail::VectorImpl +{ + public: + typedef detail::VectorImpl base_type; + + public: + VectorView(base_type & _vector, size_t n, size_t offset = 0, size_t stride = 1) + { + m_size = n; + if (stride == 1) + { + m_vector_view + = gsl_vector_subvector(_vector.get_gsl_vector(), offset, n); + m_vector = &(m_vector_view.vector); + } + else + { + m_vector_view + = gsl_vector_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n); + m_vector = &(m_vector_view.vector); + } + } + + VectorView(double* _vector, size_t n, size_t offset = 0, size_t stride = 1) + { + m_size = n; + if (stride == 1) + { + m_vector_view + = gsl_vector_view_array(_vector + offset, n); + m_vector = &(m_vector_view.vector); + } + else + { + m_vector_view + = gsl_vector_view_array_with_stride(_vector + offset, stride, n); + m_vector = &(m_vector_view.vector); + } + + } + + VectorView(const VectorView & _vector) + : base_type() + { + m_size = _vector.size(); + m_vector_view = _vector.m_vector_view; + m_vector = &(m_vector_view.vector); + } + + VectorView(Vector & _vector) + { + m_size = _vector.size(); + m_vector_view = gsl_vector_subvector(_vector.get_gsl_vector(), 0, size()); + m_vector = &(m_vector_view.vector); + } + + explicit + VectorView(gsl_vector_view _gsl_vector_view) + : m_vector_view(_gsl_vector_view) + { + m_vector = &(m_vector_view.vector); + m_size = m_vector->size; + } + + explicit + VectorView(std::vector<double> & _vector) + { + m_size = _vector.size(); + m_vector_view = gsl_vector_view_array(&(_vector[0]), _vector.size()); + m_vector = &(m_vector_view.vector); + } + + VectorView & operator=(VectorView const& _vector) + { + assert( size() == _vector.size() ); + gsl_vector_memcpy(m_vector, _vector.get_gsl_vector()); + return (*this); + } + + VectorView & operator=(base_type::base_type const& _vector) + { + assert( size() == _vector.size() ); + gsl_vector_memcpy(m_vector, _vector.get_gsl_vector()); + return (*this); + } + + VectorView & scale(double x) + { + return static_cast<VectorView&>( base_type::scale(x) ); + } + + VectorView & translate(double x) + { + return static_cast<VectorView&>( base_type::translate(x) ); + } + + VectorView & operator+=(base_type::base_type const& _vector) + { + return static_cast<VectorView&>( base_type::operator+=(_vector) ); + } + + VectorView & operator-=(base_type::base_type const& _vector) + { + return static_cast<VectorView&>( base_type::operator-=(_vector) ); + } + + friend + void swap_view(VectorView & v1, VectorView & v2); + + private: + gsl_vector_view m_vector_view; + +}; // end class VectorView + + +inline +void swap_view(VectorView & v1, VectorView & v2) +{ + assert( v1.size() == v2.size() ); + using std::swap; + swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too +} + +inline +const VectorView & const_vector_view_cast (const ConstVectorView & view) +{ + const detail::BaseVectorImpl & bvi + = static_cast<const detail::BaseVectorImpl &>(view); + const VectorView & vv = reinterpret_cast<const VectorView &>(bvi); + return vv; +} + +inline +VectorView & const_vector_view_cast (ConstVectorView & view) +{ + detail::BaseVectorImpl & bvi + = static_cast<detail::BaseVectorImpl &>(view); + VectorView & vv = reinterpret_cast<VectorView &>(bvi); + return vv; +} + + +} } // end namespaces + + +#endif /*_NL_VECTOR_H_*/ + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : |