This commit is contained in:
Jitse Niesen
2010-05-02 21:51:27 +01:00
36 changed files with 1194 additions and 883 deletions

View File

@@ -54,6 +54,8 @@ template<typename Derived> class ArrayBase
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** The base class for a given storage type. */
typedef ArrayBase StorageBaseType;
typedef ArrayBase Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl;
using ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar,
typename NumTraits<typename ei_traits<Derived>::Scalar>::Real>::operator*;

View File

@@ -2,6 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2010 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@@ -25,31 +26,55 @@
#ifndef EIGEN_GLOBAL_FUNCTIONS_H
#define EIGEN_GLOBAL_FUNCTIONS_H
#define EIGEN_ARRAY_DECLARARE_GLOBAL_UNARY(NAME,FUNCTOR) \
#define EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(NAME,FUNCTOR) \
template<typename Derived> \
inline const Eigen::CwiseUnaryOp<Eigen::FUNCTOR<typename Derived::Scalar>, Derived> \
NAME(const Eigen::ArrayBase<Derived>& x) { \
return x.derived(); \
}
#define EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(NAME,FUNCTOR) \
\
template<typename Derived> \
struct NAME##_retval<ArrayBase<Derived> > \
{ \
typedef const Eigen::CwiseUnaryOp<Eigen::FUNCTOR<typename Derived::Scalar>, Derived> type; \
}; \
template<typename Derived> \
struct NAME##_impl<ArrayBase<Derived> > \
{ \
static inline typename NAME##_retval<ArrayBase<Derived> >::type run(const Eigen::ArrayBase<Derived>& x) \
{ \
return x.derived(); \
} \
};
namespace std
{
EIGEN_ARRAY_DECLARARE_GLOBAL_UNARY(sin,ei_scalar_sin_op)
EIGEN_ARRAY_DECLARARE_GLOBAL_UNARY(cos,ei_scalar_cos_op)
EIGEN_ARRAY_DECLARARE_GLOBAL_UNARY(exp,ei_scalar_exp_op)
EIGEN_ARRAY_DECLARARE_GLOBAL_UNARY(log,ei_scalar_log_op)
EIGEN_ARRAY_DECLARARE_GLOBAL_UNARY(abs,ei_scalar_abs_op)
EIGEN_ARRAY_DECLARARE_GLOBAL_UNARY(sqrt,ei_scalar_sqrt_op)
EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(real,ei_scalar_real_op)
EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(imag,ei_scalar_imag_op)
EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(sin,ei_scalar_sin_op)
EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(cos,ei_scalar_cos_op)
EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(exp,ei_scalar_exp_op)
EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(log,ei_scalar_log_op)
EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(abs,ei_scalar_abs_op)
EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(sqrt,ei_scalar_sqrt_op)
}
namespace Eigen
{
EIGEN_ARRAY_DECLARARE_GLOBAL_UNARY(ei_sin,ei_scalar_sin_op)
EIGEN_ARRAY_DECLARARE_GLOBAL_UNARY(ei_cos,ei_scalar_cos_op)
EIGEN_ARRAY_DECLARARE_GLOBAL_UNARY(ei_exp,ei_scalar_exp_op)
EIGEN_ARRAY_DECLARARE_GLOBAL_UNARY(ei_log,ei_scalar_log_op)
EIGEN_ARRAY_DECLARARE_GLOBAL_UNARY(ei_abs,ei_scalar_abs_op)
EIGEN_ARRAY_DECLARARE_GLOBAL_UNARY(ei_sqrt,ei_scalar_sqrt_op)
EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(ei_real,ei_scalar_real_op)
EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(ei_imag,ei_scalar_imag_op)
EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(ei_sin,ei_scalar_sin_op)
EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(ei_cos,ei_scalar_cos_op)
EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(ei_exp,ei_scalar_exp_op)
EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(ei_log,ei_scalar_log_op)
EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(ei_abs,ei_scalar_abs_op)
EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(ei_abs2,ei_scalar_abs2_op)
EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(ei_sqrt,ei_scalar_sqrt_op)
}
// TODO: cleanly disable those functions that are not suppored on Array (ei_real_ref, ei_random, ei_isApprox...)
#endif // EIGEN_GLOBAL_FUNCTIONS_H

View File

@@ -5,7 +5,6 @@ ADD_SUBDIRECTORY(SVD)
ADD_SUBDIRECTORY(Cholesky)
ADD_SUBDIRECTORY(Array)
ADD_SUBDIRECTORY(Geometry)
ADD_SUBDIRECTORY(LeastSquares)
ADD_SUBDIRECTORY(Sparse)
ADD_SUBDIRECTORY(Jacobi)
ADD_SUBDIRECTORY(Householder)

View File

@@ -161,7 +161,7 @@ bool MatrixBase<Derived>::isUnitary(RealScalar prec) const
typename Derived::Nested nested(derived());
for(int i = 0; i < cols(); ++i)
{
if(!ei_isApprox(nested.col(i).squaredNorm(), static_cast<Scalar>(1), prec))
if(!ei_isApprox(nested.col(i).squaredNorm(), static_cast<RealScalar>(1), prec))
return false;
for(int j = 0; j < i; ++j)
if(!ei_isMuchSmallerThan(nested.col(i).dot(nested.col(j)), static_cast<Scalar>(1), prec))

View File

@@ -180,7 +180,7 @@ struct ei_functor_traits<ei_scalar_quotient_op<Scalar> > {
Cost = 2 * NumTraits<Scalar>::MulCost,
PacketAccess = ei_packet_traits<Scalar>::size>1
#if (defined EIGEN_VECTORIZE)
&& NumTraits<Scalar>::HasFloatingPoint
&& !NumTraits<Scalar>::IsInteger
#endif
};
};
@@ -384,7 +384,7 @@ template<typename Scalar1,typename Scalar2>
struct ei_functor_traits<ei_scalar_multiple2_op<Scalar1,Scalar2> >
{ enum { Cost = NumTraits<Scalar1>::MulCost, PacketAccess = false }; };
template<typename Scalar, bool HasFloatingPoint>
template<typename Scalar, bool IsInteger>
struct ei_scalar_quotient1_impl {
typedef typename ei_packet_traits<Scalar>::type PacketScalar;
// FIXME default copy constructors seems bugged with std::complex<>
@@ -396,11 +396,11 @@ struct ei_scalar_quotient1_impl {
const Scalar m_other;
};
template<typename Scalar>
struct ei_functor_traits<ei_scalar_quotient1_impl<Scalar,true> >
struct ei_functor_traits<ei_scalar_quotient1_impl<Scalar,false> >
{ enum { Cost = NumTraits<Scalar>::MulCost, PacketAccess = ei_packet_traits<Scalar>::size>1 }; };
template<typename Scalar>
struct ei_scalar_quotient1_impl<Scalar,false> {
struct ei_scalar_quotient1_impl<Scalar,true> {
// FIXME default copy constructors seems bugged with std::complex<>
EIGEN_STRONG_INLINE ei_scalar_quotient1_impl(const ei_scalar_quotient1_impl& other) : m_other(other.m_other) { }
EIGEN_STRONG_INLINE ei_scalar_quotient1_impl(const Scalar& other) : m_other(other) {}
@@ -408,7 +408,7 @@ struct ei_scalar_quotient1_impl<Scalar,false> {
typename ei_makeconst<typename NumTraits<Scalar>::Nested>::type m_other;
};
template<typename Scalar>
struct ei_functor_traits<ei_scalar_quotient1_impl<Scalar,false> >
struct ei_functor_traits<ei_scalar_quotient1_impl<Scalar,true> >
{ enum { Cost = 2 * NumTraits<Scalar>::MulCost, PacketAccess = false }; };
/** \internal
@@ -420,13 +420,13 @@ struct ei_functor_traits<ei_scalar_quotient1_impl<Scalar,false> >
* \sa class CwiseUnaryOp, MatrixBase::operator/
*/
template<typename Scalar>
struct ei_scalar_quotient1_op : ei_scalar_quotient1_impl<Scalar, NumTraits<Scalar>::HasFloatingPoint > {
struct ei_scalar_quotient1_op : ei_scalar_quotient1_impl<Scalar, NumTraits<Scalar>::IsInteger > {
EIGEN_STRONG_INLINE ei_scalar_quotient1_op(const Scalar& other)
: ei_scalar_quotient1_impl<Scalar, NumTraits<Scalar>::HasFloatingPoint >(other) {}
: ei_scalar_quotient1_impl<Scalar, NumTraits<Scalar>::IsInteger >(other) {}
};
template<typename Scalar>
struct ei_functor_traits<ei_scalar_quotient1_op<Scalar> >
: ei_functor_traits<ei_scalar_quotient1_impl<Scalar, NumTraits<Scalar>::HasFloatingPoint> >
: ei_functor_traits<ei_scalar_quotient1_impl<Scalar, NumTraits<Scalar>::IsInteger> >
{};
// nullary functors

View File

@@ -26,6 +26,8 @@
#ifndef EIGEN_FUZZY_H
#define EIGEN_FUZZY_H
// TODO support small integer types properly i.e. do exact compare on coeffs --- taking a HS norm is guaranteed to cause integer overflow.
#ifndef EIGEN_LEGACY_COMPARES
/** \returns \c true if \c *this is approximately equal to \a other, within the precision

View File

@@ -126,8 +126,8 @@ DenseBase<Derived>::format(const IOFormat& fmt) const
return WithFormat<Derived>(derived(), fmt);
}
template<typename Scalar>
struct ei_significant_decimals_impl
template<typename Scalar, bool IsInteger>
struct ei_significant_decimals_default_impl
{
typedef typename NumTraits<Scalar>::Real RealScalar;
static inline int run()
@@ -136,6 +136,20 @@ struct ei_significant_decimals_impl
}
};
template<typename Scalar>
struct ei_significant_decimals_default_impl<Scalar, true>
{
static inline int run()
{
return 0;
}
};
template<typename Scalar>
struct ei_significant_decimals_impl
: ei_significant_decimals_default_impl<Scalar, NumTraits<Scalar>::IsInteger>
{};
/** \internal
* print the matrix \a _m to the output stream \a s using the output format \a fmt */
template<typename Derived>
@@ -153,13 +167,13 @@ std::ostream & ei_print_matrix(std::ostream & s, const Derived& _m, const IOForm
}
else if(fmt.precision == FullPrecision)
{
if (NumTraits<Scalar>::HasFloatingPoint)
if (NumTraits<Scalar>::IsInteger)
{
explicit_precision = ei_significant_decimals_impl<Scalar>::run();
explicit_precision = 0;
}
else
{
explicit_precision = 0;
explicit_precision = ei_significant_decimals_impl<Scalar>::run();
}
}
else

File diff suppressed because it is too large Load Diff

View File

@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2006-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@@ -27,157 +27,121 @@
/** \class NumTraits
*
* \brief Holds some data about the various numeric (i.e. scalar) types allowed by Eigen.
* \brief Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
*
* \param T the numeric type about which this class provides data. Recall that Eigen allows
* only the following types for \a T: \c int, \c float, \c double,
* \c std::complex<float>, \c std::complex<double>, and \c long \c double (especially
* useful to enforce x87 arithmetics when SSE is the default).
* \param T the numeric type at hand
*
* The provided data consists of everything that is supported by std::numeric_limits, plus:
* This class stores enums, typedefs and static methods giving information about a numeric type.
*
* The provided data consists of:
* \li A typedef \a Real, giving the "real part" type of \a T. If \a T is already real,
* then \a Real is just a typedef to \a T. If \a T is \c std::complex<U> then \a Real
* is a typedef to \a U.
* \li A typedef \a FloatingPoint, giving the "floating-point type" of \a T. If \a T is
* \c int, then \a FloatingPoint is a typedef to \c double. Otherwise, \a FloatingPoint
* is a typedef to \a T.
* \li A typedef \a NonInteger, giving the type that should be used for operations producing non-integral values,
* such as quotients, square roots, etc. If \a T is a floating-point type, then this typedef just gives
* \a T again. Note however that many Eigen functions such as ei_sqrt simply refuse to
* take integers. Outside of a few cases, Eigen doesn't do automatic type promotion. Thus, this typedef is
* only intended as a helper for code that needs to explicitly promote types.
* \li A typedef \a Nested giving the type to use to nest a value inside of the expression tree. If you don't know what
* this means, just use \a T here.
* \li An enum value \a IsComplex. It is equal to 1 if \a T is a \c std::complex
* type, and to 0 otherwise.
* \li An enum \a HasFloatingPoint. It is equal to \c 0 if \a T is \c int,
* and to \c 1 otherwise.
* \li An enum value \a IsInteger. It is equal to \c 1 if \a T is an integer type such as \c int,
* and to \c 0 otherwise.
* \li Enum values ReadCost, AddCost and MulCost representing a rough estimate of the number of CPU cycles needed
* to by move / add / mul instructions respectively, assuming the data is already stored in CPU registers.
* Stay vague here. No need to do architecture-specific stuff.
* \li An enum value \a IsSigned. It is equal to \c 1 if \a T is a signed type and to 0 if \a T is unsigned.
* \li An epsilon() function which, unlike std::numeric_limits::epsilon(), returns a \a Real instead of a \a T.
* \li A dummy_precision() function returning a weak epsilon value. It is mainly used by the fuzzy comparison operators.
* \li Two highest() and lowest() functions returning the highest and lowest possible values respectively.
* \li A dummy_precision() function returning a weak epsilon value. It is mainly used as a default
* value by the fuzzy comparison operators.
* \li highest() and lowest() functions returning the highest and lowest possible values respectively.
*/
template<typename T> struct NumTraits;
template<typename T> struct ei_default_float_numtraits
: std::numeric_limits<T>
template<typename T> struct GenericNumTraits
{
inline static T highest() { return std::numeric_limits<T>::max(); }
inline static T lowest() { return -std::numeric_limits<T>::max(); }
};
enum {
IsInteger = std::numeric_limits<T>::is_integer,
IsSigned = std::numeric_limits<T>::is_signed,
IsComplex = 0,
ReadCost = 1,
AddCost = 1,
MulCost = 1
};
template<typename T> struct ei_default_integral_numtraits
: std::numeric_limits<T>
{
inline static T dummy_precision() { return T(0); }
typedef T Real;
typedef typename ei_meta_if<
IsInteger,
typename ei_meta_if<sizeof(T)<=2, float, double>::ret,
T
>::ret NonInteger;
typedef T Nested;
inline static Real epsilon() { return std::numeric_limits<T>::epsilon(); }
inline static Real dummy_precision()
{
// make sure to override this for floating-point types
return Real(0);
}
inline static T highest() { return std::numeric_limits<T>::max(); }
inline static T lowest() { return std::numeric_limits<T>::min(); }
};
template<> struct NumTraits<int>
: ei_default_integral_numtraits<int>
{
typedef int Real;
typedef double FloatingPoint;
typedef int Nested;
enum {
IsComplex = 0,
HasFloatingPoint = 0,
ReadCost = 1,
AddCost = 1,
MulCost = 1
};
};
template<typename T> struct NumTraits : GenericNumTraits<T>
{};
template<> struct NumTraits<float>
: ei_default_float_numtraits<float>
: GenericNumTraits<float>
{
typedef float Real;
typedef float FloatingPoint;
typedef float Nested;
enum {
IsComplex = 0,
HasFloatingPoint = 1,
ReadCost = 1,
AddCost = 1,
MulCost = 1
};
inline static float dummy_precision() { return 1e-5f; }
};
template<> struct NumTraits<double>
: ei_default_float_numtraits<double>
template<> struct NumTraits<double> : GenericNumTraits<double>
{
typedef double Real;
typedef double FloatingPoint;
typedef double Nested;
enum {
IsComplex = 0,
HasFloatingPoint = 1,
ReadCost = 1,
AddCost = 1,
MulCost = 1
};
inline static double dummy_precision() { return 1e-12; }
};
template<> struct NumTraits<long double>
: GenericNumTraits<long double>
{
static inline long double dummy_precision() { return 1e-15l; }
};
template<typename _Real> struct NumTraits<std::complex<_Real> >
: ei_default_float_numtraits<std::complex<_Real> >
: GenericNumTraits<std::complex<_Real> >
{
typedef _Real Real;
typedef std::complex<_Real> FloatingPoint;
typedef std::complex<_Real> Nested;
enum {
IsComplex = 1,
HasFloatingPoint = NumTraits<Real>::HasFloatingPoint,
ReadCost = 2,
AddCost = 2 * NumTraits<Real>::AddCost,
MulCost = 4 * NumTraits<Real>::MulCost + 2 * NumTraits<Real>::AddCost
};
inline static Real epsilon() { return std::numeric_limits<Real>::epsilon(); }
inline static Real epsilon() { return NumTraits<Real>::epsilon(); }
inline static Real dummy_precision() { return NumTraits<Real>::dummy_precision(); }
};
template<> struct NumTraits<long long int>
: ei_default_integral_numtraits<long long int>
template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
struct NumTraits<Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols> >
{
typedef long long int Real;
typedef long double FloatingPoint;
typedef long long int Nested;
typedef Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols> ArrayType;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Array<RealScalar, Rows, Cols, Options, MaxRows, MaxCols> Real;
typedef typename NumTraits<Scalar>::NonInteger NonIntegerScalar;
typedef Array<NonIntegerScalar, Rows, Cols, Options, MaxRows, MaxCols> NonInteger;
typedef ArrayType & Nested;
enum {
IsComplex = 0,
HasFloatingPoint = 0,
ReadCost = 1,
AddCost = 1,
MulCost = 1
IsComplex = NumTraits<Scalar>::IsComplex,
IsInteger = NumTraits<Scalar>::IsInteger,
IsSigned = NumTraits<Scalar>::IsSigned,
ReadCost = ArrayType::SizeAtCompileTime==Dynamic ? Dynamic : ArrayType::SizeAtCompileTime * NumTraits<Scalar>::ReadCost,
AddCost = ArrayType::SizeAtCompileTime==Dynamic ? Dynamic : ArrayType::SizeAtCompileTime * NumTraits<Scalar>::AddCost,
MulCost = ArrayType::SizeAtCompileTime==Dynamic ? Dynamic : ArrayType::SizeAtCompileTime * NumTraits<Scalar>::MulCost
};
};
template<> struct NumTraits<long double>
: ei_default_float_numtraits<long double>
{
typedef long double Real;
typedef long double FloatingPoint;
typedef long double Nested;
enum {
IsComplex = 0,
HasFloatingPoint = 1,
ReadCost = 1,
AddCost = 1,
MulCost = 1
};
static inline long double dummy_precision() { return NumTraits<double>::dummy_precision(); }
};
template<> struct NumTraits<bool>
: ei_default_integral_numtraits<bool>
{
typedef bool Real;
typedef float FloatingPoint;
typedef bool Nested;
enum {
IsComplex = 0,
HasFloatingPoint = 0,
ReadCost = 1,
AddCost = 1,
MulCost = 1
};
};
#endif // EIGEN_NUMTRAITS_H

View File

@@ -133,9 +133,11 @@ inline Derived& DenseBase<Derived>::operator*=(const Scalar& other)
template<typename Derived>
inline Derived& DenseBase<Derived>::operator/=(const Scalar& other)
{
SelfCwiseBinaryOp<typename ei_meta_if<NumTraits<Scalar>::HasFloatingPoint,ei_scalar_product_op<Scalar>,ei_scalar_quotient_op<Scalar> >::ret, Derived> tmp(derived());
SelfCwiseBinaryOp<typename ei_meta_if<NumTraits<Scalar>::IsInteger,
ei_scalar_quotient_op<Scalar>,
ei_scalar_product_op<Scalar> >::ret, Derived> tmp(derived());
typedef typename Derived::PlainObject PlainObject;
tmp = PlainObject::Constant(rows(),cols(), NumTraits<Scalar>::HasFloatingPoint ? Scalar(1)/other : other);
tmp = PlainObject::Constant(rows(),cols(), NumTraits<Scalar>::IsInteger ? other : Scalar(1)/other);
return derived();
}

View File

@@ -65,7 +65,7 @@
YOU_CALLED_A_FIXED_SIZE_METHOD_ON_A_DYNAMIC_SIZE_MATRIX_OR_VECTOR,
YOU_CALLED_A_DYNAMIC_SIZE_METHOD_ON_A_FIXED_SIZE_MATRIX_OR_VECTOR,
UNALIGNED_LOAD_AND_STORE_OPERATIONS_UNIMPLEMENTED_ON_ALTIVEC,
NUMERIC_TYPE_MUST_BE_FLOATING_POINT,
THIS_FUNCTION_IS_NOT_FOR_INTEGER_NUMERIC_TYPES,
NUMERIC_TYPE_MUST_BE_REAL,
COEFFICIENT_WRITE_ACCESS_TO_SELFADJOINT_NOT_SUPPORTED,
WRITING_TO_TRIANGULAR_PART_WITH_UNIT_DIAGONAL_IS_NOT_SUPPORTED,
@@ -158,6 +158,9 @@
) \
)
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE) \
EIGEN_STATIC_ASSERT(!NumTraits<TYPE>::IsInteger, THIS_FUNCTION_IS_NOT_FOR_INTEGER_NUMERIC_TYPES)
// static assertion failing if it is guaranteed at compile-time that the two matrix expression types have different sizes
#define EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(TYPE0,TYPE1) \
EIGEN_STATIC_ASSERT( \

View File

@@ -46,7 +46,7 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim)
typedef _Scalar Scalar;
typedef NumTraits<Scalar> ScalarTraits;
typedef typename ScalarTraits::Real RealScalar;
typedef typename ScalarTraits::FloatingPoint FloatingPoint;
typedef typename ScalarTraits::NonInteger NonInteger;
typedef Matrix<Scalar,AmbientDimAtCompileTime,1> VectorType;
/** Define constants to name the corners of a 1D, 2D or 3D axis aligned bounding box */
@@ -174,11 +174,10 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim)
VectorType r;
for(int d=0; d<dim(); ++d)
{
if(ScalarTraits::HasFloatingPoint)
if(!ScalarTraits::IsInteger)
{
r[d] = m_min[d] + (m_max[d]-m_min[d])
* (ei_random<Scalar>() + ei_random_amplitude<Scalar>())
/ (Scalar(2)*ei_random_amplitude<Scalar>() );
* ei_random<Scalar>(Scalar(0), Scalar(1));
}
else
r[d] = ei_random(m_min[d], m_max[d]);
@@ -260,15 +259,15 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim)
* \sa squaredExteriorDistance()
*/
template<typename Derived>
inline FloatingPoint exteriorDistance(const MatrixBase<Derived>& p) const
{ return ei_sqrt(FloatingPoint(squaredExteriorDistance(p))); }
inline NonInteger exteriorDistance(const MatrixBase<Derived>& p) const
{ return ei_sqrt(NonInteger(squaredExteriorDistance(p))); }
/** \returns the distance between the boxes \a b and \c *this,
* and zero if the boxes intersect.
* \sa squaredExteriorDistance()
*/
inline FloatingPoint exteriorDistance(const AlignedBox& b) const
{ return ei_sqrt(FloatingPoint(squaredExteriorDistance(b))); }
inline NonInteger exteriorDistance(const AlignedBox& b) const
{ return ei_sqrt(NonInteger(squaredExteriorDistance(b))); }
/** \returns \c *this with scalar type casted to \a NewScalarType
*

View File

@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2008-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@@ -325,7 +325,7 @@ struct ei_inverse_impl : public ReturnByValue<ei_inverse_impl<MatrixType> >
template<typename Derived>
inline const ei_inverse_impl<Derived> MatrixBase<Derived>::inverse() const
{
EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT)
EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsInteger,THIS_FUNCTION_IS_NOT_FOR_INTEGER_NUMERIC_TYPES)
ei_assert(rows() == cols());
return ei_inverse_impl<Derived>(derived());
}

View File

@@ -1,6 +0,0 @@
FILE(GLOB Eigen_LeastSquares_SRCS "*.h")
INSTALL(FILES
${Eigen_LeastSquares_SRCS}
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/LeastSquares COMPONENT Devel
)

View File

@@ -1,181 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_LEASTSQUARES_H
#define EIGEN_LEASTSQUARES_H
/** \ingroup LeastSquares_Module
*
* \leastsquares_module
*
* For a set of points, this function tries to express
* one of the coords as a linear (affine) function of the other coords.
*
* This is best explained by an example. This function works in full
* generality, for points in a space of arbitrary dimension, and also over
* the complex numbers, but for this example we will work in dimension 3
* over the real numbers (doubles).
*
* So let us work with the following set of 5 points given by their
* \f$(x,y,z)\f$ coordinates:
* @code
Vector3d points[5];
points[0] = Vector3d( 3.02, 6.89, -4.32 );
points[1] = Vector3d( 2.01, 5.39, -3.79 );
points[2] = Vector3d( 2.41, 6.01, -4.01 );
points[3] = Vector3d( 2.09, 5.55, -3.86 );
points[4] = Vector3d( 2.58, 6.32, -4.10 );
* @endcode
* Suppose that we want to express the second coordinate (\f$y\f$) as a linear
* expression in \f$x\f$ and \f$z\f$, that is,
* \f[ y=ax+bz+c \f]
* for some constants \f$a,b,c\f$. Thus, we want to find the best possible
* constants \f$a,b,c\f$ so that the plane of equation \f$y=ax+bz+c\f$ fits
* best the five above points. To do that, call this function as follows:
* @code
Vector3d coeffs; // will store the coefficients a, b, c
linearRegression(
5,
&points,
&coeffs,
1 // the coord to express as a function of
// the other ones. 0 means x, 1 means y, 2 means z.
);
* @endcode
* Now the vector \a coeffs is approximately
* \f$( 0.495 , -1.927 , -2.906 )\f$.
* Thus, we get \f$a=0.495, b = -1.927, c = -2.906\f$. Let us check for
* instance how near points[0] is from the plane of equation \f$y=ax+bz+c\f$.
* Looking at the coords of points[0], we see that:
* \f[ax+bz+c = 0.495 * 3.02 + (-1.927) * (-4.32) + (-2.906) = 6.91.\f]
* On the other hand, we have \f$y=6.89\f$. We see that the values
* \f$6.91\f$ and \f$6.89\f$
* are near, so points[0] is very near the plane of equation \f$y=ax+bz+c\f$.
*
* Let's now describe precisely the parameters:
* @param numPoints the number of points
* @param points the array of pointers to the points on which to perform the linear regression
* @param result pointer to the vector in which to store the result.
This vector must be of the same type and size as the
data points. The meaning of its coords is as follows.
For brevity, let \f$n=Size\f$,
\f$r_i=result[i]\f$,
and \f$f=funcOfOthers\f$. Denote by
\f$x_0,\ldots,x_{n-1}\f$
the n coordinates in the n-dimensional space.
Then the resulting equation is:
\f[ x_f = r_0 x_0 + \cdots + r_{f-1}x_{f-1}
+ r_{f+1}x_{f+1} + \cdots + r_{n-1}x_{n-1} + r_n. \f]
* @param funcOfOthers Determines which coord to express as a function of the
others. Coords are numbered starting from 0, so that a
value of 0 means \f$x\f$, 1 means \f$y\f$,
2 means \f$z\f$, ...
*
* \sa fitHyperplane()
*/
template<typename VectorType>
void linearRegression(int numPoints,
VectorType **points,
VectorType *result,
int funcOfOthers )
{
typedef typename VectorType::Scalar Scalar;
typedef Hyperplane<Scalar, VectorType::SizeAtCompileTime> HyperplaneType;
const int size = points[0]->size();
result->resize(size);
HyperplaneType h(size);
fitHyperplane(numPoints, points, &h);
for(int i = 0; i < funcOfOthers; i++)
result->coeffRef(i) = - h.coeffs()[i] / h.coeffs()[funcOfOthers];
for(int i = funcOfOthers; i < size; i++)
result->coeffRef(i) = - h.coeffs()[i+1] / h.coeffs()[funcOfOthers];
}
/** \ingroup LeastSquares_Module
*
* \leastsquares_module
*
* This function is quite similar to linearRegression(), so we refer to the
* documentation of this function and only list here the differences.
*
* The main difference from linearRegression() is that this function doesn't
* take a \a funcOfOthers argument. Instead, it finds a general equation
* of the form
* \f[ r_0 x_0 + \cdots + r_{n-1}x_{n-1} + r_n = 0, \f]
* where \f$n=Size\f$, \f$r_i=retCoefficients[i]\f$, and we denote by
* \f$x_0,\ldots,x_{n-1}\f$ the n coordinates in the n-dimensional space.
*
* Thus, the vector \a retCoefficients has size \f$n+1\f$, which is another
* difference from linearRegression().
*
* In practice, this function performs an hyper-plane fit in a total least square sense
* via the following steps:
* 1 - center the data to the mean
* 2 - compute the covariance matrix
* 3 - pick the eigenvector corresponding to the smallest eigenvalue of the covariance matrix
* The ratio of the smallest eigenvalue and the second one gives us a hint about the relevance
* of the solution. This value is optionally returned in \a soundness.
*
* \sa linearRegression()
*/
template<typename VectorType, typename HyperplaneType>
void fitHyperplane(int numPoints,
VectorType **points,
HyperplaneType *result,
typename NumTraits<typename VectorType::Scalar>::Real* soundness = 0)
{
typedef typename VectorType::Scalar Scalar;
typedef Matrix<Scalar,VectorType::SizeAtCompileTime,VectorType::SizeAtCompileTime> CovMatrixType;
EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType)
ei_assert(numPoints >= 1);
int size = points[0]->size();
ei_assert(size+1 == result->coeffs().size());
// compute the mean of the data
VectorType mean = VectorType::Zero(size);
for(int i = 0; i < numPoints; ++i)
mean += *(points[i]);
mean /= numPoints;
// compute the covariance matrix
CovMatrixType covMat = CovMatrixType::Zero(size, size);
for(int i = 0; i < numPoints; ++i)
{
VectorType diff = (*(points[i]) - mean).conjugate();
covMat += diff * diff.adjoint();
}
// now we just have to pick the eigen vector with smallest eigen value
SelfAdjointEigenSolver<CovMatrixType> eig(covMat);
result->normal() = eig.eigenvectors().col(0);
if (soundness)
*soundness = eig.eigenvalues().coeff(0)/eig.eigenvalues().coeff(1);
// let's compute the constant coefficient such that the
// plane pass trough the mean point:
result->offset() = - (result->normal().cwiseProduct(mean)).sum();
}
#endif // EIGEN_LEASTSQUARES_H