diff options
author | Benoit Jacob <jacob.benoit.1@gmail.com> | 2007-10-10 06:09:56 +0000 |
---|---|---|
committer | Benoit Jacob <jacob.benoit.1@gmail.com> | 2007-10-10 06:09:56 +0000 |
commit | 06e1e0d83be56b4b38041fb19749ae79303c09c5 (patch) | |
tree | c41f29aeae71c11f5102bcf14e785b46046df011 /src/Core | |
parent | 7f0a546a817e97e0fc55ae21b7db811315dd94fe (diff) |
fix dot product, add norm/norm2/normalized
add fuzzy compares for matrices/vectors
add random matrix/vector generation
Diffstat (limited to 'src/Core')
-rw-r--r-- | src/Core/Dot.h | 31 | ||||
-rw-r--r-- | src/Core/Fuzzy.h | 81 | ||||
-rw-r--r-- | src/Core/Numeric.h | 66 | ||||
-rw-r--r-- | src/Core/Object.h | 14 | ||||
-rw-r--r-- | src/Core/Random.h | 67 | ||||
-rw-r--r-- | src/Core/ScalarOps.h | 2 | ||||
-rw-r--r-- | src/Core/Util.h | 1 |
7 files changed, 217 insertions, 45 deletions
diff --git a/src/Core/Dot.h b/src/Core/Dot.h index 175c3e57d..0d000ea70 100644 --- a/src/Core/Dot.h +++ b/src/Core/Dot.h @@ -31,23 +31,17 @@ struct EiDotUnroller { static void run(const Derived1 &v1, const Derived2& v2, typename Derived1::Scalar &dot) { - const int i = Index - 1; EiDotUnroller<Index-1, Size, Derived1, Derived2>::run(v1, v2, dot); - if(i == Size - 1) - dot = v1[i] * EiConj(v2[i]); - else - dot += v1[i] * EiConj(v2[i]); + dot += v1[Index-1] * EiConj(v2[Index-1]); } }; template<int Size, typename Derived1, typename Derived2> -struct EiDotUnroller<0, Size, Derived1, Derived2> +struct EiDotUnroller<1, Size, Derived1, Derived2> { static void run(const Derived1 &v1, const Derived2& v2, typename Derived1::Scalar &dot) { - EI_UNUSED(v1); - EI_UNUSED(v2); - EI_UNUSED(dot); + dot = v1[0] * EiConj(v2[0]); } }; @@ -80,4 +74,23 @@ Scalar EiObject<Scalar, Derived>::dot(const OtherDerived& other) const return res; } +template<typename Scalar, typename Derived> +typename EiNumTraits<Scalar>::Real EiObject<Scalar, Derived>::norm2() const +{ + assert(IsVector); + return EiReal(dot(*this)); +} + +template<typename Scalar, typename Derived> +typename EiNumTraits<Scalar>::Real EiObject<Scalar, Derived>::norm() const +{ + return EiSqrt(norm2()); +} + +template<typename Scalar, typename Derived> +EiScalarProduct<Derived> EiObject<Scalar, Derived>::normalized() const +{ + return (*this) / norm(); +} + #endif // EI_DOT_H diff --git a/src/Core/Fuzzy.h b/src/Core/Fuzzy.h new file mode 100644 index 000000000..44b20615d --- /dev/null +++ b/src/Core/Fuzzy.h @@ -0,0 +1,81 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2006-2007 Benoit Jacob <jacob@math.jussieu.fr> +// +// Eigen is free software; 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 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 General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along +// with Eigen; if not, write to the Free Software Foundation, Inc., 51 +// Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. +// +// As a special exception, if other files instantiate templates or use macros +// or functions from this file, or you compile this file and link it +// with other works to produce a work based on this file, this file does not +// by itself cause the resulting work to be covered by the GNU General Public +// License. This exception does not invalidate any other reasons why a work +// based on this file might be covered by the GNU General Public License. + +#ifndef EI_FUZZY_H +#define EI_FUZZY_H + +template<typename Scalar, typename Derived> +template<typename OtherDerived> +bool EiObject<Scalar, Derived>::isApprox(const OtherDerived& other) const +{ + if(IsVector) + { + return((*this - other).norm2() + <= std::min(norm2(), other.norm2()) + * EiAbs2(EiNumTraits<Scalar>::epsilon())); + } + else + { + for(int i = 0; i < cols(); i++) + if(!col(i).isApprox(other.col(i))) + return false; + return true; + } +} + +template<typename Scalar, typename Derived> +bool EiObject<Scalar, Derived>::isNegligble(const Scalar& other) const +{ + if(IsVector) + { + return(norm2() <= EiAbs2(other) * EiAbs2(EiNumTraits<Scalar>::epsilon())); + } + else + { + for(int i = 0; i < cols(); i++) + if(!col(i).isNegligible(other)) + return false; + return true; + } +} + +template<typename Scalar, typename Derived> +template<typename OtherDerived> +bool EiObject<Scalar, Derived>::isNegligble(const EiObject<Scalar, OtherDerived>& other) const +{ + if(IsVector) + { + return(norm2() <= other.norm2() * EiAbs2(EiNumTraits<Scalar>::epsilon())); + } + else + { + for(int i = 0; i < cols(); i++) + if(!col(i).isNegligible(other.col(i))) + return false; + return true; + } +} + +#endif // EI_FUZZY_H
\ No newline at end of file diff --git a/src/Core/Numeric.h b/src/Core/Numeric.h index 48a36f236..23a6c03e6 100644 --- a/src/Core/Numeric.h +++ b/src/Core/Numeric.h @@ -27,9 +27,9 @@ #ifndef EI_NUMERIC_H #define EI_NUMERIC_H -template<typename T> struct EiTraits; +template<typename T> struct EiNumTraits; -template<> struct EiTraits<int> +template<> struct EiNumTraits<int> { typedef int Real; typedef double FloatingPoint; @@ -45,16 +45,16 @@ template<> struct EiTraits<int> static double sqrt(const int& x) { return std::sqrt(static_cast<double>(x)); } static int abs(const int& x) { return std::abs(x); } static int abs2(const int& x) { return x*x; } - static int random() + static int rand() { // "rand()%21" would be bad. always use the high-order bits, not the low-order bits. // note: here (gcc 4.1) static_cast<int> seems to round the nearest int. // I don't know if that's part of the standard. - return -10 + static_cast<int>(rand() / ((RAND_MAX + 1.0)/20.0)); + return -10 + static_cast<int>(std::rand() / ((RAND_MAX + 1.0)/20.0)); } }; -template<> struct EiTraits<float> +template<> struct EiNumTraits<float> { typedef float Real; typedef float FloatingPoint; @@ -70,13 +70,13 @@ template<> struct EiTraits<float> static float sqrt(const float& x) { return std::sqrt(x); } static float abs(const float& x) { return std::abs(x); } static float abs2(const float& x) { return x*x; } - static float random() + static float rand() { - return rand() / (RAND_MAX/20.0f) - 10.0f; + return std::rand() / (RAND_MAX/20.0f) - 10.0f; } }; -template<> struct EiTraits<double> +template<> struct EiNumTraits<double> { typedef double Real; typedef double FloatingPoint; @@ -92,23 +92,23 @@ template<> struct EiTraits<double> static double sqrt(const double& x) { return std::sqrt(x); } static double abs(const double& x) { return std::abs(x); } static double abs2(const double& x) { return x*x; } - static double random() + static double rand() { - return rand() / (RAND_MAX/20.0) - 10.0; + return std::rand() / (RAND_MAX/20.0) - 10.0; } }; -template<typename _Real> struct EiTraits<std::complex<_Real> > +template<typename _Real> struct EiNumTraits<std::complex<_Real> > { typedef _Real Real; typedef std::complex<Real> Complex; typedef std::complex<double> FloatingPoint; - typedef typename EiTraits<Real>::FloatingPoint RealFloatingPoint; + typedef typename EiNumTraits<Real>::FloatingPoint RealFloatingPoint; static const bool IsComplex = true; - static const bool HasFloatingPoint = EiTraits<Real>::HasFloatingPoint; + static const bool HasFloatingPoint = EiNumTraits<Real>::HasFloatingPoint; - static Real epsilon() { return EiTraits<Real>::epsilon(); } + static Real epsilon() { return EiNumTraits<Real>::epsilon(); } static Real real(const Complex& x) { return std::real(x); } static Real imag(const Complex& x) { return std::imag(x); } static Complex conj(const Complex& x) { return std::conj(x); } @@ -118,49 +118,49 @@ template<typename _Real> struct EiTraits<std::complex<_Real> > { return std::abs(static_cast<FloatingPoint>(x)); } static Real abs2(const Complex& x) { return std::real(x) * std::real(x) + std::imag(x) * std::imag(x); } - static Complex random() + static Complex rand() { - return Complex(EiTraits<Real>::random(), EiTraits<Real>::random()); + return Complex(EiNumTraits<Real>::rand(), EiNumTraits<Real>::rand()); } }; -template<typename T> typename EiTraits<T>::Real EiReal(const T& x) -{ return EiTraits<T>::real(x); } +template<typename T> typename EiNumTraits<T>::Real EiReal(const T& x) +{ return EiNumTraits<T>::real(x); } -template<typename T> typename EiTraits<T>::Real EiImag(const T& x) -{ return EiTraits<T>::imag(x); } +template<typename T> typename EiNumTraits<T>::Real EiImag(const T& x) +{ return EiNumTraits<T>::imag(x); } template<typename T> T EiConj(const T& x) -{ return EiTraits<T>::conj(x); } +{ return EiNumTraits<T>::conj(x); } -template<typename T> typename EiTraits<T>::FloatingPoint EiSqrt(const T& x) -{ return EiTraits<T>::sqrt(x); } +template<typename T> typename EiNumTraits<T>::FloatingPoint EiSqrt(const T& x) +{ return EiNumTraits<T>::sqrt(x); } -template<typename T> typename EiTraits<T>::RealFloatingPoint EiAbs(const T& x) -{ return EiTraits<T>::abs(x); } +template<typename T> typename EiNumTraits<T>::RealFloatingPoint EiAbs(const T& x) +{ return EiNumTraits<T>::abs(x); } -template<typename T> typename EiTraits<T>::Real EiAbs2(const T& x) -{ return EiTraits<T>::abs2(x); } +template<typename T> typename EiNumTraits<T>::Real EiAbs2(const T& x) +{ return EiNumTraits<T>::abs2(x); } -template<typename T> T EiRandom() -{ return EiTraits<T>::random(); } +template<typename T> T EiRand() +{ return EiNumTraits<T>::rand(); } template<typename T> bool EiNegligible(const T& a, const T& b) { - return(EiAbs(a) <= EiAbs(b) * EiTraits<T>::epsilon()); + return(EiAbs(a) <= EiAbs(b) * EiNumTraits<T>::epsilon()); } template<typename T> bool EiApprox(const T& a, const T& b) { - if(EiTraits<T>::IsFloat) - return(EiAbs(a - b) <= std::min(EiAbs(a), EiAbs(b)) * EiTraits<T>::epsilon()); + if(EiNumTraits<T>::IsFloat) + return(EiAbs(a - b) <= std::min(EiAbs(a), EiAbs(b)) * EiNumTraits<T>::epsilon()); else return(a == b); } template<typename T> bool EiLessThanOrApprox(const T& a, const T& b) { - if(EiTraits<T>::IsFloat) + if(EiNumTraits<T>::IsFloat) return(a < b || EiApprox(a, b)); else return(a <= b); diff --git a/src/Core/Object.h b/src/Core/Object.h index 89c62b7b5..d1a310a20 100644 --- a/src/Core/Object.h +++ b/src/Core/Object.h @@ -42,6 +42,7 @@ template<typename Scalar, typename Derived> class EiObject typedef typename EiForwardDecl<Derived>::Ref Ref; typedef typename EiForwardDecl<Derived>::ConstRef ConstRef; + typedef typename EiNumTraits<Scalar>::Real RealScalar; int rows() const { return static_cast<const Derived *>(this)->_rows(); } int cols() const { return static_cast<const Derived *>(this)->_cols(); } @@ -92,8 +93,17 @@ template<typename Scalar, typename Derived> class EiObject template<typename OtherDerived> Scalar dot(const OtherDerived& other) const; - Scalar norm2() const { assert(IsVector); return dot(*this); } - Scalar norm() const { assert(IsVector); return EiSqrt(dot(*this)); } + RealScalar norm2() const; + RealScalar norm() const; + EiScalarProduct<Derived> normalized() const; + + static EiEval<EiRandom<Derived> > + random(int rows = RowsAtCompileTime, int cols = ColsAtCompileTime); + + template<typename OtherDerived> + bool isApprox(const OtherDerived& other) const; + template<typename OtherDerived> + bool isNegligible(const OtherDerived& other) const; template<typename OtherDerived> EiMatrixProduct<Derived, OtherDerived> diff --git a/src/Core/Random.h b/src/Core/Random.h new file mode 100644 index 000000000..70485028f --- /dev/null +++ b/src/Core/Random.h @@ -0,0 +1,67 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2006-2007 Benoit Jacob <jacob@math.jussieu.fr> +// +// Eigen is free software; 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 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 General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along +// with Eigen; if not, write to the Free Software Foundation, Inc., 51 +// Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. +// +// As a special exception, if other files instantiate templates or use macros +// or functions from this file, or you compile this file and link it +// with other works to produce a work based on this file, this file does not +// by itself cause the resulting work to be covered by the GNU General Public +// License. This exception does not invalidate any other reasons why a work +// based on this file might be covered by the GNU General Public License. + +#ifndef EI_RANDOM_H +#define EI_RANDOM_H + +template<typename MatrixType> class EiRandom + : public EiObject<typename MatrixType::Scalar, EiRandom<MatrixType> > +{ + public: + typedef typename MatrixType::Scalar Scalar; + friend class EiObject<Scalar, EiRandom<MatrixType> >; + + static const int RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime; + + EiRandom(int rows, int cols) : m_rows(rows), m_cols(cols) + { + assert(rows > 0 && cols > 0); + } + + private: + EiRandom& _ref() { return *this; } + const EiRandom& _constRef() const { return *this; } + int _rows() const { return m_rows; } + int _cols() const { return m_cols; } + + Scalar _read(int row, int col) const + { + EI_UNUSED(row); + EI_UNUSED(col); + return EiRand<Scalar>(); + } + + protected: + int m_rows, m_cols; +}; + +template<typename Scalar, typename Derived> +EiEval<EiRandom<Derived> > EiObject<Scalar, Derived>::random(int rows, int cols) +{ + return EiRandom<Derived>(rows, cols).eval(); +} + +#endif // EI_RANDOM_H diff --git a/src/Core/ScalarOps.h b/src/Core/ScalarOps.h index 2f56ff3b4..c058978e6 100644 --- a/src/Core/ScalarOps.h +++ b/src/Core/ScalarOps.h @@ -83,7 +83,7 @@ EiScalarProduct<Derived> \ operator/(const EiObject<Scalar, Derived>& matrix, \ OtherScalar scalar) \ { \ - assert(EiTraits<Scalar>::HasFloatingPoint); \ + assert(EiNumTraits<Scalar>::HasFloatingPoint); \ return matrix * (static_cast<Scalar>(1) / scalar); \ } \ \ diff --git a/src/Core/Util.h b/src/Core/Util.h index 859913908..460eecb1a 100644 --- a/src/Core/Util.h +++ b/src/Core/Util.h @@ -54,6 +54,7 @@ template<typename Lhs, typename Rhs> class EiSum; template<typename Lhs, typename Rhs> class EiDifference; template<typename Lhs, typename Rhs> class EiMatrixProduct; template<typename MatrixType> class EiScalarProduct; +template<typename MatrixType> class EiRandom; template<typename ExpressionType> class EiEval; template<typename T> struct EiForwardDecl |