|
[Rivet-svn] r1986 - in trunk: . include include/Rivet include/Rivet/Math include/Rivet/Math/eigen src src/Testblackhole at projects.hepforge.org blackhole at projects.hepforge.orgTue Nov 3 14:34:59 GMT 2009
Author: buckley Date: Tue Nov 3 14:34:58 2009 New Revision: 1986 Log: Removing eigen and replacing with eigen2: expression-template-based linear algebra library which will hopefully give a bit of a performance boost to matrix operations. I've added more matrix and vector tests to the unit tests, and they all pass, but please keep an eye out for regressions in observable plots. Deleted: trunk/include/Rivet/Math/eigen/ Modified: trunk/ChangeLog trunk/configure.ac trunk/include/Makefile.am trunk/include/Rivet/Makefile.am trunk/include/Rivet/Math/MathHeader.hh trunk/include/Rivet/Math/MathUtils.hh trunk/include/Rivet/Math/Matrix3.hh trunk/include/Rivet/Math/MatrixN.hh trunk/include/Rivet/Math/StdHeader.hh trunk/include/Rivet/Math/VectorN.hh trunk/src/Run.cc trunk/src/Test/testMatVec.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Tue Nov 3 14:29:40 2009 (r1985) +++ trunk/ChangeLog Tue Nov 3 14:34:58 2009 (r1986) @@ -1,3 +1,11 @@ +2009-11-03 Andy Buckley <andy at insectnation.org> + + * Adding more assertion checks to linear algebra testing. + + * Replacing eigen linear algebra library with eigen2, which is + based on expression templates and may give some performance + benefits. + 2009-11-02 Hendrik Hoeth <hendrik.hoeth at cern.ch> * Fixing normalisation issue with stacked histograms in make-plots. Modified: trunk/configure.ac ============================================================================== --- trunk/configure.ac Tue Nov 3 14:29:40 2009 (r1985) +++ trunk/configure.ac Tue Nov 3 14:34:58 2009 (r1986) @@ -230,6 +230,7 @@ ## Set default build flags AM_CPPFLAGS="-I\$(top_srcdir)/include -I\$(top_builddir)/include" +AM_CPPFLAGS="$AM_CPPFLAGS -I\$(top_srcdir)/include/eigen2" AM_CPPFLAGS="$AM_CPPFLAGS -I\$(BOOSTINCPATH) \$(BOOST_CPPFLAGS)" AM_CPPFLAGS="$AM_CPPFLAGS -I\$(GSLINCPATH)" AM_CPPFLAGS="$AM_CPPFLAGS -I\$(HEPMCINCPATH)" @@ -237,9 +238,7 @@ AC_CEDAR_CHECKCXXFLAG([-pedantic], [AM_CXXFLAGS="$AM_CXXFLAGS -pedantic "]) AC_CEDAR_CHECKCXXFLAG([-ansi], [AM_CXXFLAGS="$AM_CXXFLAGS -ansi "]) AC_CEDAR_CHECKCXXFLAG([-Wall], [AM_CXXFLAGS="$AM_CXXFLAGS -Wall "]) -dnl if test x$with_root == xyes; then -dnl AC_CEDAR_CHECKCXXFLAG([-Wno-long-long], [AM_CXXFLAGS="$AM_CXXFLAGS -Wno-long-long "]) -dnl fi +AC_CEDAR_CHECKCXXFLAG([-Wno-long-long], [AM_CXXFLAGS="$AM_CXXFLAGS -Wno-long-long "]) ## Debug flag (default=none) Modified: trunk/include/Makefile.am ============================================================================== --- trunk/include/Makefile.am Tue Nov 3 14:29:40 2009 (r1985) +++ trunk/include/Makefile.am Tue Nov 3 14:34:58 2009 (r1986) @@ -1,7 +1,9 @@ SUBDIRS = Rivet +EXTRA_DIST = eigen2 + ## Remove these when YODA is available... -EXTRA_DIST = TinyXML +EXTRA_DIST += TinyXML nobase_include_HEADERS = \ LWH/AIAnalysisFactory.h \ LWH/AITreeFactory.h \ Modified: trunk/include/Rivet/Makefile.am ============================================================================== --- trunk/include/Rivet/Makefile.am Tue Nov 3 14:29:40 2009 (r1985) +++ trunk/include/Rivet/Makefile.am Tue Nov 3 14:34:58 2009 (r1986) @@ -107,16 +107,4 @@ Math/Vector4.hh \ Math/Math.hh \ Math/Units.hh \ - Math/Constants.hh \ - Math/eigen/util.h \ - Math/eigen/regressioninternal.h \ - Math/eigen/regression.h \ - Math/eigen/vector.h \ - Math/eigen/ludecompositionbase.h \ - Math/eigen/ludecomposition.h \ - Math/eigen/linearsolver.h \ - Math/eigen/linearsolverbase.h \ - Math/eigen/matrix.h \ - Math/eigen/vectorbase.h \ - Math/eigen/projective.h \ - Math/eigen/matrixbase.h + Math/Constants.hh Modified: trunk/include/Rivet/Math/MathHeader.hh ============================================================================== --- trunk/include/Rivet/Math/MathHeader.hh Tue Nov 3 14:29:40 2009 (r1985) +++ trunk/include/Rivet/Math/MathHeader.hh Tue Nov 3 14:34:58 2009 (r1986) @@ -1,5 +1,5 @@ -#ifndef RIVET_Math_MathHeader -#define RIVET_Math_MathHeader +#ifndef RIVET_MATH_MATHHEADER +#define RIVET_MATH_MATHHEADER #include <stdexcept> #include <string> @@ -10,9 +10,10 @@ #include <map> #include <vector> #include <algorithm> +#include <cassert> -#include "Rivet/Math/eigen/vector.h" -#include "Rivet/Math/eigen/matrix.h" +// #include "Rivet/Math/eigen/vector.h" +// #include "Rivet/Math/eigen/matrix.h" namespace Rivet { Modified: trunk/include/Rivet/Math/MathUtils.hh ============================================================================== --- trunk/include/Rivet/Math/MathUtils.hh Tue Nov 3 14:29:40 2009 (r1985) +++ trunk/include/Rivet/Math/MathUtils.hh Tue Nov 3 14:34:58 2009 (r1986) @@ -3,6 +3,7 @@ #define RIVET_MathUtils_HH #include "Rivet/Math/MathHeader.hh" +#include "Rivet/RivetBoost.hh" namespace Rivet { Modified: trunk/include/Rivet/Math/Matrix3.hh ============================================================================== --- trunk/include/Rivet/Math/Matrix3.hh Tue Nov 3 14:29:40 2009 (r1985) +++ trunk/include/Rivet/Math/Matrix3.hh Tue Nov 3 14:34:58 2009 (r1986) @@ -6,6 +6,8 @@ #include "Rivet/Math/MatrixN.hh" #include "Rivet/Math/Vector3.hh" +#include <Eigen/Geometry> + namespace Rivet { @@ -16,8 +18,13 @@ Matrix3(const Matrix<3>& m3) : Matrix<3>::Matrix<3>(m3) { } Matrix3(const Vector3& axis, const double angle) { - const Vector3 normaxis = axis.unit(); - _matrix.loadRotation3(angle, normaxis._vec); + if (Rivet::isZero(angle)) { + _matrix.setIdentity(); + } else { + const Vector3 normaxis = axis.unit(); + Eigen::AngleAxis<double> eaxis(angle, normaxis._vec); + _matrix = Eigen::Transform<double,3>(eaxis).linear(); + } } Matrix3(const Vector3& from, const Vector3& to) { @@ -41,10 +48,11 @@ Matrix3& setAsRotation(const Vector3& from, const Vector3& to) { const double theta = angle(from, to); if (Rivet::isZero(theta)) { - _matrix.loadIdentity(); + _matrix.setIdentity(); } else { const Vector3 normaxis = cross(from, to).unit(); - _matrix.loadRotation3(theta, normaxis._vec); + Eigen::AngleAxis<double> eaxis(theta, normaxis._vec); + _matrix = Eigen::Transform<double,3>(eaxis).linear(); } return *this; } Modified: trunk/include/Rivet/Math/MatrixN.hh ============================================================================== --- trunk/include/Rivet/Math/MatrixN.hh Tue Nov 3 14:29:40 2009 (r1985) +++ trunk/include/Rivet/Math/MatrixN.hh Tue Nov 3 14:34:58 2009 (r1986) @@ -5,6 +5,9 @@ #include "Rivet/Math/MathUtils.hh" #include "Rivet/Math/Vectors.hh" +#include <Eigen/Core> +#include <Eigen/LU> + namespace Rivet { @@ -62,7 +65,7 @@ public: Matrix() { - _matrix.loadZero(); + _matrix.setZero(); } Matrix(const Matrix<N>& other) { @@ -102,7 +105,7 @@ } Vector<N> getColumn(const size_t col) const { - const Eigen::Vector<double,N> eVec = _matrix.column(col); + const EVector eVec = _matrix.col(col); Vector<N> rtn; for (size_t i = 0; i < N; ++i) { rtn.set(i, _matrix(i,col)); @@ -119,12 +122,12 @@ Matrix<N> transpose() const { Matrix<N> tmp = *this; - tmp._matrix.replaceWithAdjoint(); + tmp._matrix.transposeInPlace(); return tmp; } // Matrix<N>& transposeInPlace() { - // _matrix.replaceWithAdjoint(); + // _matrix.transposeInPlace(); // return *this; // } @@ -239,7 +242,8 @@ } protected: - typedef Eigen::Matrix<double,N> EMatrix; + typedef Eigen::Matrix<double,N,1> EVector; + typedef Eigen::Matrix<double,N,N> EMatrix; EMatrix _matrix; }; @@ -376,6 +380,20 @@ } +template <size_t N> +inline bool fuzzyEquals(const Matrix<N>& ma, const Matrix<N>& mb, double tolerance=1E-5) { + for (size_t i = 0; i < N; ++i) { + for (size_t j = 0; j < N; ++j) { + const double a = ma.get(i, j); + const double b = mb.get(i, j); + if (!Rivet::fuzzyEquals(a, b, tolerance)) return false; + } + } + return true; +} + + } + #endif Modified: trunk/include/Rivet/Math/StdHeader.hh ============================================================================== --- trunk/include/Rivet/Math/StdHeader.hh Tue Nov 3 14:29:40 2009 (r1985) +++ trunk/include/Rivet/Math/StdHeader.hh Tue Nov 3 14:34:58 2009 (r1986) @@ -10,8 +10,8 @@ #include <map> #include <vector> -#include "eigen/vector.h" -#include "eigen/matrix.h" +// #include "eigen/vector.h" +// #include "eigen/matrix.h" using std::string; using std::ostream; Modified: trunk/include/Rivet/Math/VectorN.hh ============================================================================== --- trunk/include/Rivet/Math/VectorN.hh Tue Nov 3 14:29:40 2009 (r1985) +++ trunk/include/Rivet/Math/VectorN.hh Tue Nov 3 14:34:58 2009 (r1986) @@ -4,6 +4,8 @@ #include "Rivet/Math/MathHeader.hh" #include "Rivet/Math/MathUtils.hh" +#include <Eigen/Core> + namespace Rivet { @@ -23,21 +25,22 @@ friend Vector<M> multiply(const Matrix<M>& a, const Vector<M>& b); public: - Vector() { _vec.loadZero(); } + Vector() { _vec.setZero(); } Vector(const Vector<N>& other) : _vec(other._vec) { } - const double& get(const size_t index) const { + /// Get vector elements + const double get(const size_t index) const { if (index >= N) { throw std::runtime_error("Tried to access an invalid vector index."); } else { - return _vec(index); + return _vec[index]; } } /// Direct access to vector elements by index. - const double& operator[](const size_t index) const { + const double operator[](const size_t index) const { return get(index); } @@ -122,12 +125,13 @@ if (index >= N) { throw std::runtime_error("Tried to access an invalid vector index."); } else { - return _vec(index); + return _vec[index]; } } /// Vector - Eigen::Vector<double,N> _vec; + typedef Eigen::Matrix<double,N,1> EVector; + EVector _vec; }; @@ -149,6 +153,7 @@ ///////////////////////////////////////////////// + template <size_t N> inline const string toString(const Vector<N>& v) { ostringstream out; @@ -168,6 +173,20 @@ } + ///////////////////////////////////////////////// + + + template <size_t N> + inline bool fuzzyEquals(const Vector<N>& va, const Vector<N>& vb, double tolerance=1E-5) { + for (size_t i = 0; i < N; ++i) { + const double a = va.get(i); + const double b = vb.get(i); + if (!Rivet::fuzzyEquals(a, b, tolerance)) return false; + } + return true; + } + + } #endif Modified: trunk/src/Run.cc ============================================================================== --- trunk/src/Run.cc Tue Nov 3 14:29:40 2009 (r1985) +++ trunk/src/Run.cc Tue Nov 3 14:34:58 2009 (r1986) @@ -72,7 +72,7 @@ _ah.removeIncompatibleAnalyses(beamIds(*evt)); if (num_anas_requested > 0 && _ah.analysisNames().size() == 0) { Log::getLog("Rivet.Run") << Log::ERROR - << "All analyses were incompatible with the first event's beams" + << "All analyses were incompatible with the first event's beams\n" << "Exiting, since this probably isn't intentional!" << endl; delete evt; return false; Modified: trunk/src/Test/testMatVec.cc ============================================================================== --- trunk/src/Test/testMatVec.cc Tue Nov 3 14:29:40 2009 (r1985) +++ trunk/src/Test/testMatVec.cc Tue Nov 3 14:34:58 2009 (r1986) @@ -6,24 +6,30 @@ #include "Rivet/Math/Matrices.hh" using namespace std; -using namespace Rivet; int main() { - + using namespace Rivet; + FourVector a(1,0,0,0); cout << a << ": interval = " << a.invariant() << endl; + assert(fuzzyEquals(a.invariant(), 1)); a.setZ(1); + assert(isZero(a.invariant())); cout << a << ": interval = " << a.invariant() << endl; a.setY(2).setZ(3); cout << a << ": interval = " << a.invariant() << endl; + assert(fuzzyEquals(a.invariant(), -12)); cout << a << ": vector = " << a.vector3() << endl << endl; - FourMomentum b(1,0,0,1); + FourMomentum b(1,0,0,0); cout << b << ": mass = " << b.mass() << endl; + assert(fuzzyEquals(b.mass2(), 1)); b.setPz(1); cout << b << ": mass = " << b.mass() << endl; + assert(isZero(b.mass2())); b.setPy(2).setPz(3).setE(6); cout << b << ": mass = " << b.mass2() << endl; + assert(fuzzyEquals(b.mass2(), 23)); cout << b << ": vector = " << b.vector3() << endl << endl; Matrix3 m; @@ -53,14 +59,51 @@ m3.setRow(0, Vector3(2,3,0)).setRow(1, Vector3(1,4,3)).setRow(2, Vector3(0,1,2)); cout << m1+m2 << " == " << m3 << ": " << (m1+m2 == m3 ? "true" : "false") << endl; cout << endl; + Vector3 v3(1,2,3); cout << "Vector: " << v3 << endl; cout << "Invert: " << v3 << " --> " << -v3 << endl; const Matrix3 rot90(Vector3(0,0,1), PI/2.0); - cout << "Rot 90: " << v3 << " --90deg--> " << rot90*v3 << endl; - cout << "Rot 2 x 90: " << v3 << " --90deg--> " << rot90*rot90*v3 << endl; - cout << "Rot 90*-90: "<< v3 << " --90deg--> " << rot90*rot90.inverse()*v3 << endl; + const Matrix3 rot90m(Vector3(0,0,1), -PI/2.0); + const Matrix3 rot180(Vector3(0,0,1), PI); + const Matrix3 rot180m(Vector3(0,0,1), -PI); + const Vector3 v3_90 = rot90*v3; + cout << "Rot 90: " << v3 << " ---> " << v3_90 << endl; + const Vector3 v3_90m = rot90m*v3; + cout << "Rot -90: " << v3 << " ---> " << v3_90m << endl; + const Vector3 v3_180 = rot180*v3; + cout << "Rot 180: " << v3 << " ---> " << v3_180 << endl; + const Vector3 v3_180m = rot180m*v3; + cout << "Rot -180: " << v3 << " ---> " << v3_180m << endl; + assert(fuzzyEquals(v3_180, v3_180m)); + + const Vector3 v3_9090 = rot90*rot90*v3; + cout << "Rot 2 x 90: " << v3 << " ---> " << v3_9090 << endl; + assert(fuzzyEquals(v3_180, v3_9090)); + + const Vector3 v3_90m90m = rot90m*rot90m*v3; + cout << "Rot 2 x -90: " << v3 << " ---> " << v3_90m90m << endl; + assert(fuzzyEquals(v3_180, v3_90m90m)); + + const Vector3 v3_9090m = rot90*rot90m*v3; + const Vector3 v3_90m90 = rot90m*rot90*v3; + cout << "Rot 90*-90: "<< v3 << " ---> " << v3_9090m << endl; + cout << "Rot -90*90: "<< v3 << " ---> " << v3_90m90 << endl; + assert(fuzzyEquals(v3, v3_9090m)); + assert(fuzzyEquals(v3, v3_90m90)); + + const Vector3 v3_90i = rot90.inverse()*v3; + cout << "Rot (90)^-1: "<< v3 << " ---> " << v3_90i << endl; + assert(fuzzyEquals(v3_90i, v3_90m)); + + const Vector3 v3_9090i = rot90*rot90.inverse()*v3; + const Vector3 v3_90i90 = rot90.inverse()*rot90*v3; + cout << "Rot 90*(90)^-1: "<< v3 << " ---> " << v3_9090i << endl; + cout << "Rot (90)^-1*90: "<< v3 << " ---> " << v3_90i90 << endl; + assert(fuzzyEquals(v3, v3_9090i)); + assert(fuzzyEquals(v3, v3_90i90)); + const Matrix3 rot1(Vector3(0,1,0), PI/180.0); cout << "Rot 0 x 45 x 1: " << v3 << endl; for (size_t i = 0; i < 8; ++i) { @@ -69,7 +112,7 @@ } cout << "Rot " << i+1 << " x 45 x 1: " << v3 << endl; } - assert(v3 == Vector3(1,2,3)); + assert(fuzzyEquals(v3, Vector3(1,2,3))); cout << endl; cout << "Boosts:" << endl;
More information about the Rivet-svn mailing list |