|
|
[Rivet-svn] r1992 - in trunk: . include include/Rivet/Math include/Rivet/Math/eigen include/eigen2 pyext
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Wed Nov 4 00:46:05 GMT 2009
Author: buckley
Date: Wed Nov 4 00:46:05 2009
New Revision: 1992
Log:
Reverting to use of bundled eigen1, until upgrade to eigen2 can be done without the dependency being a pain
Added:
trunk/include/Rivet/Math/eigen/
- copied from r1985, trunk/include/Rivet/Math/eigen/
Deleted:
trunk/include/eigen2/
Modified:
trunk/ChangeLog
trunk/include/Makefile.am
trunk/include/Rivet/Math/MathHeader.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/pyext/setup.py.in
Modified: trunk/ChangeLog
==============================================================================
--- trunk/ChangeLog Wed Nov 4 00:12:46 2009 (r1991)
+++ trunk/ChangeLog Wed Nov 4 00:46:05 2009 (r1992)
@@ -2,10 +2,6 @@
* 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/include/Makefile.am
==============================================================================
--- trunk/include/Makefile.am Wed Nov 4 00:12:46 2009 (r1991)
+++ trunk/include/Makefile.am Wed Nov 4 00:46:05 2009 (r1992)
@@ -1,9 +1,7 @@
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 \
@@ -31,3 +29,5 @@
LWH/Tree.h \
LWH/Axis.h \
LWH/VariAxis.h
+
+#EXTRA_DIST += eigen2
Modified: trunk/include/Rivet/Math/MathHeader.hh
==============================================================================
--- trunk/include/Rivet/Math/MathHeader.hh Wed Nov 4 00:12:46 2009 (r1991)
+++ trunk/include/Rivet/Math/MathHeader.hh Wed Nov 4 00:46:05 2009 (r1992)
@@ -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,10 +10,6 @@
#include <map>
#include <vector>
#include <algorithm>
-#include <cassert>
-
-// #include "Rivet/Math/eigen/vector.h"
-// #include "Rivet/Math/eigen/matrix.h"
namespace Rivet {
Modified: trunk/include/Rivet/Math/Matrix3.hh
==============================================================================
--- trunk/include/Rivet/Math/Matrix3.hh Wed Nov 4 00:12:46 2009 (r1991)
+++ trunk/include/Rivet/Math/Matrix3.hh Wed Nov 4 00:46:05 2009 (r1992)
@@ -6,8 +6,6 @@
#include "Rivet/Math/MatrixN.hh"
#include "Rivet/Math/Vector3.hh"
-#include <Eigen/Geometry>
-
namespace Rivet {
@@ -18,13 +16,8 @@
Matrix3(const Matrix<3>& m3) : Matrix<3>::Matrix<3>(m3) { }
Matrix3(const Vector3& axis, const double angle) {
- 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();
- }
+ const Vector3 normaxis = axis.unit();
+ _matrix.loadRotation3(angle, normaxis._vec);
}
Matrix3(const Vector3& from, const Vector3& to) {
@@ -48,11 +41,10 @@
Matrix3& setAsRotation(const Vector3& from, const Vector3& to) {
const double theta = angle(from, to);
if (Rivet::isZero(theta)) {
- _matrix.setIdentity();
+ _matrix.loadIdentity();
} else {
const Vector3 normaxis = cross(from, to).unit();
- Eigen::AngleAxis<double> eaxis(theta, normaxis._vec);
- _matrix = Eigen::Transform<double,3>(eaxis).linear();
+ _matrix.loadRotation3(theta, normaxis._vec);
}
return *this;
}
Modified: trunk/include/Rivet/Math/MatrixN.hh
==============================================================================
--- trunk/include/Rivet/Math/MatrixN.hh Wed Nov 4 00:12:46 2009 (r1991)
+++ trunk/include/Rivet/Math/MatrixN.hh Wed Nov 4 00:46:05 2009 (r1992)
@@ -5,412 +5,414 @@
#include "Rivet/Math/MathUtils.hh"
#include "Rivet/Math/Vectors.hh"
-#include <Eigen/Core>
-#include <Eigen/LU>
+#include "Rivet/Math/eigen/matrix.h"
namespace Rivet {
-template <size_t N>
-class Matrix;
-typedef Matrix<4> Matrix4;
-
-template <size_t N>
-Matrix<N> multiply(const Matrix<N>& a, const Matrix<N>& b);
-template <size_t N>
-Matrix<N> divide(const Matrix<N>&, const double);
-template <size_t N>
-Matrix<N> operator*(const Matrix<N>& a, const Matrix<N>& b);
-
-
-///////////////////////////////////
-
-
-template <size_t N>
-class Matrix {
- template <size_t M>
- friend Matrix<M> add(const Matrix<M>&, const Matrix<M>&);
- template <size_t M>
- friend Matrix<M> multiply(const double, const Matrix<M>&);
- template <size_t M>
- friend Matrix<M> multiply(const Matrix<M>&, const Matrix<M>&);
- template <size_t M>
- friend Vector<M> multiply(const Matrix<M>&, const Vector<M>&);
- template <size_t M>
- friend Matrix<M> divide(const Matrix<M>&, const double);
-
-public:
- static Matrix<N> mkZero() {
- Matrix<N> rtn;
- return rtn;
- }
-
- static Matrix<N> mkDiag(Vector<N> diag) {
- Matrix<N> rtn;
- for (size_t i = 0; i < N; ++i) {
- rtn.set(i, i, diag[i]);
+ template <size_t N>
+ class Matrix;
+ typedef Matrix<4> Matrix4;
+
+ template <size_t N>
+ Matrix<N> multiply(const Matrix<N>& a, const Matrix<N>& b);
+ template <size_t N>
+ Matrix<N> divide(const Matrix<N>&, const double);
+ template <size_t N>
+ Matrix<N> operator*(const Matrix<N>& a, const Matrix<N>& b);
+
+
+ ///////////////////////////////////
+
+
+ template <size_t N>
+ class Matrix {
+ template <size_t M>
+ friend Matrix<M> add(const Matrix<M>&, const Matrix<M>&);
+ template <size_t M>
+ friend Matrix<M> multiply(const double, const Matrix<M>&);
+ template <size_t M>
+ friend Matrix<M> multiply(const Matrix<M>&, const Matrix<M>&);
+ template <size_t M>
+ friend Vector<M> multiply(const Matrix<M>&, const Vector<M>&);
+ template <size_t M>
+ friend Matrix<M> divide(const Matrix<M>&, const double);
+
+ public:
+ static Matrix<N> mkZero() {
+ Matrix<N> rtn;
+ return rtn;
+ }
+
+ static Matrix<N> mkDiag(Vector<N> diag) {
+ Matrix<N> rtn;
+ for (size_t i = 0; i < N; ++i) {
+ rtn.set(i, i, diag[i]);
+ }
+ return rtn;
}
- return rtn;
- }
-
- static Matrix<N> mkIdentity() {
- Matrix<N> rtn;
- for (size_t i = 0; i < N; ++i) {
- rtn.set(i, i, 1);
+
+ static Matrix<N> mkIdentity() {
+ Matrix<N> rtn;
+ for (size_t i = 0; i < N; ++i) {
+ rtn.set(i, i, 1);
+ }
+ return rtn;
}
- return rtn;
- }
-public:
+ public:
- Matrix() {
- _matrix.setZero();
- }
+ Matrix() {
+ _matrix.loadZero();
+ }
- Matrix(const Matrix<N>& other) {
- _matrix = other._matrix;
- }
-
- Matrix& set(const size_t i, const size_t j, const double value) {
- if (i < N && j < N) {
- _matrix(i, j) = value;
- } else {
- throw std::runtime_error("Attempted set access outside matrix bounds.");
+ Matrix(const Matrix<N>& other) {
+ _matrix = other._matrix;
}
- return *this;
- }
-
- double get(const size_t i, const size_t j) const {
- if (i < N && j < N) {
- return _matrix(i, j);
- } else {
- throw std::runtime_error("Attempted get access outside matrix bounds.");
+
+ Matrix& set(const size_t i, const size_t j, const double value) {
+ if (i < N && j < N) {
+ _matrix(i, j) = value;
+ } else {
+ throw std::runtime_error("Attempted set access outside matrix bounds.");
+ }
+ return *this;
}
- }
- Vector<N> getRow(const size_t row) const {
- Vector<N> rtn;
- for (size_t i = 0; i < N; ++i) {
- rtn.set(i, _matrix(row,i));
+ double get(const size_t i, const size_t j) const {
+ if (i < N && j < N) {
+ return _matrix(i, j);
+ } else {
+ throw std::runtime_error("Attempted get access outside matrix bounds.");
+ }
}
- return rtn;
- }
- Matrix<N>& setRow(const size_t row, const Vector<N>& r) {
- for (size_t i = 0; i < N; ++i) {
- _matrix(row,i) = r.get(i);
+ Vector<N> getRow(const size_t row) const {
+ Vector<N> rtn;
+ for (size_t i = 0; i < N; ++i) {
+ rtn.set(i, _matrix(row,i));
+ }
+ return rtn;
}
- return *this;
- }
- Vector<N> getColumn(const size_t col) const {
- const EVector eVec = _matrix.col(col);
- Vector<N> rtn;
- for (size_t i = 0; i < N; ++i) {
- rtn.set(i, _matrix(i,col));
+ Matrix<N>& setRow(const size_t row, const Vector<N>& r) {
+ for (size_t i = 0; i < N; ++i) {
+ _matrix(row,i) = r.get(i);
+ }
+ return *this;
}
- return rtn;
- }
- Matrix<N>& setColumn(const size_t col, const Vector<N>& c) {
- for (size_t i = 0; i < N; ++i) {
- _matrix(i,col) = c.get(i);
+ Vector<N> getColumn(const size_t col) const {
+ const Eigen::Vector<double,N> eVec = _matrix.column(col);
+ Vector<N> rtn;
+ for (size_t i = 0; i < N; ++i) {
+ rtn.set(i, _matrix(i,col));
+ }
+ return rtn;
}
- return *this;
- }
- Matrix<N> transpose() const {
- Matrix<N> tmp = *this;
- tmp._matrix.transposeInPlace();
- return tmp;
- }
+ Matrix<N>& setColumn(const size_t col, const Vector<N>& c) {
+ for (size_t i = 0; i < N; ++i) {
+ _matrix(i,col) = c.get(i);
+ }
+ return *this;
+ }
- // Matrix<N>& transposeInPlace() {
- // _matrix.transposeInPlace();
- // return *this;
- // }
+ Matrix<N> transpose() const {
+ Matrix<N> tmp = *this;
+ tmp._matrix.replaceWithAdjoint();
+ return tmp;
+ }
- /// Calculate inverse
- Matrix<N> inverse() const {
- Matrix<N> tmp;
- tmp._matrix = _matrix.inverse();
- return tmp;
- }
+ // Matrix<N>& transposeInPlace() {
+ // _matrix.replaceWithAdjoint();
+ // return *this;
+ // }
- /// Calculate determinant
- double det() const {
- return _matrix.determinant();
- }
+ /// Calculate inverse
+ Matrix<N> inverse() const {
+ Matrix<N> tmp;
+ tmp._matrix = _matrix.inverse();
+ return tmp;
+ }
- /// Calculate trace
- double trace() const {
- return _matrix.trace();
- }
+ /// Calculate determinant
+ double det() const {
+ return _matrix.determinant();
+ }
- /// Negate
- Matrix<N> operator-() const {
- Matrix<N> rtn;
- rtn._matrix = -_matrix;
- return rtn;
- }
+ /// Calculate trace
+ double trace() const {
+ double tr = 0.0;
+ for (size_t i = 0; i < N; ++i) {
+ tr += _matrix(i,i);
+ }
+ return tr;
+ // return _matrix.trace();
+ }
- /// Get dimensionality
- const size_t size() const {
- return N;
- }
+ /// Negate
+ Matrix<N> operator-() const {
+ Matrix<N> rtn;
+ rtn._matrix = -_matrix;
+ return rtn;
+ }
+
+ /// Get dimensionality
+ const size_t size() const {
+ return N;
+ }
- /// Index-wise check for nullness, allowing for numerical precision
- const bool isZero(double tolerance=1E-5) const {
- for (size_t i=0; i < N; ++i) {
- for (size_t j=0; j < N; ++j) {
- if (! Rivet::isZero(_matrix(i,j), tolerance) ) return false;
+ /// Index-wise check for nullness, allowing for numerical precision
+ const bool isZero(double tolerance=1E-5) const {
+ for (size_t i=0; i < N; ++i) {
+ for (size_t j=0; j < N; ++j) {
+ if (! Rivet::isZero(_matrix(i,j), tolerance) ) return false;
+ }
}
+ return true;
}
- return true;
- }
- /// Check for index-wise equality, allowing for numerical precision
- const bool isEqual(Matrix<N> other) const {
- for (size_t i=0; i < N; ++i) {
- for (size_t j=i; j < N; ++j) {
- if (! Rivet::fuzzyEquals(_matrix(i,j), other._matrix(i,j)) ) return false;
+ /// Check for index-wise equality, allowing for numerical precision
+ const bool isEqual(Matrix<N> other) const {
+ for (size_t i=0; i < N; ++i) {
+ for (size_t j=i; j < N; ++j) {
+ if (! Rivet::isZero(_matrix(i,j) - other._matrix(i,j)) ) return false;
+ }
}
+ return true;
}
- return true;
- }
- /// Check for symmetry under transposition
- const bool isSymm() const {
- return isEqual(this->transpose());
- }
+ /// Check for symmetry under transposition
+ const bool isSymm() const {
+ return isEqual(this->transpose());
+ }
- /// Check that all off-diagonal elements are zero, allowing for numerical precision
- const bool isDiag() const {
- for (size_t i=0; i < N; ++i) {
- for (size_t j=0; j < N; ++j) {
- if (i == j) continue;
- if (! Rivet::isZero(_matrix(i,j)) ) return false;
+ /// Check that all off-diagonal elements are zero, allowing for numerical precision
+ const bool isDiag() const {
+ for (size_t i=0; i < N; ++i) {
+ for (size_t j=0; j < N; ++j) {
+ if (i == j) continue;
+ if (! Rivet::isZero(_matrix(i,j)) ) return false;
+ }
}
+ return true;
}
- return true;
- }
- bool operator==(const Matrix<N>& a) const {
- return _matrix == a._matrix;
- }
+ bool operator==(const Matrix<N>& a) const {
+ return _matrix == a._matrix;
+ }
- bool operator!=(const Matrix<N>& a) const {
- return _matrix != a._matrix;
- }
+ bool operator!=(const Matrix<N>& a) const {
+ return _matrix != a._matrix;
+ }
- bool operator<(const Matrix<N>& a) const {
- return _matrix < a._matrix;
- }
+ bool operator<(const Matrix<N>& a) const {
+ return _matrix < a._matrix;
+ }
- bool operator<=(const Matrix<N>& a) const {
- return _matrix <= a._matrix;
- }
+ bool operator<=(const Matrix<N>& a) const {
+ return _matrix <= a._matrix;
+ }
- bool operator>(const Matrix<N>& a) const {
- return _matrix > a._matrix;
- }
+ bool operator>(const Matrix<N>& a) const {
+ return _matrix > a._matrix;
+ }
- bool operator>=(const Matrix<N>& a) const {
- return _matrix >= a._matrix;
- }
+ bool operator>=(const Matrix<N>& a) const {
+ return _matrix >= a._matrix;
+ }
- Matrix<N>& operator*=(const Matrix<N>& m) {
- _matrix = _matrix * m._matrix;
- return *this;
- }
+ Matrix<N>& operator*=(const Matrix<N>& m) {
+ _matrix = _matrix * m._matrix;
+ return *this;
+ }
- Matrix<N>& operator*=(const double a) {
- _matrix *= a;
- return *this;
- }
+ Matrix<N>& operator*=(const double a) {
+ _matrix *= a;
+ return *this;
+ }
- Matrix<N>& operator/=(const double a) {
- _matrix /= a;
- return *this;
- }
+ Matrix<N>& operator/=(const double a) {
+ _matrix /= a;
+ return *this;
+ }
- Matrix<N>& operator+=(const Matrix<N>& m) {
- _matrix += m._matrix;
- return *this;
- }
+ Matrix<N>& operator+=(const Matrix<N>& m) {
+ _matrix += m._matrix;
+ return *this;
+ }
- Matrix<N>& operator-=(const Matrix<N>& m) {
- _matrix -= m._matrix;
- return *this;
- }
+ Matrix<N>& operator-=(const Matrix<N>& m) {
+ _matrix -= m._matrix;
+ return *this;
+ }
-protected:
- typedef Eigen::Matrix<double,N,1> EVector;
- typedef Eigen::Matrix<double,N,N> EMatrix;
- EMatrix _matrix;
-};
+ protected:
+ typedef Eigen::Matrix<double,N> EMatrix;
+ EMatrix _matrix;
+ };
-/////////////////////////////////
+ /////////////////////////////////
-template <size_t N>
-inline Matrix<N> add(const Matrix<N>& a, const Matrix<N>& b) {
- Matrix<N> result;
- result._matrix = a._matrix + b._matrix;
- return result;
-}
+ template <size_t N>
+ inline Matrix<N> add(const Matrix<N>& a, const Matrix<N>& b) {
+ Matrix<N> result;
+ result._matrix = a._matrix + b._matrix;
+ return result;
+ }
-template <size_t N>
-inline Matrix<N> subtract(const Matrix<N>& a, const Matrix<N>& b) {
- return add(a, -b);
-}
+ template <size_t N>
+ inline Matrix<N> subtract(const Matrix<N>& a, const Matrix<N>& b) {
+ return add(a, -b);
+ }
-template <size_t N>
-inline Matrix<N> operator+(const Matrix<N> a, const Matrix<N>& b) {
- return add(a, b);
-}
+ template <size_t N>
+ inline Matrix<N> operator+(const Matrix<N> a, const Matrix<N>& b) {
+ return add(a, b);
+ }
-template <size_t N>
-inline Matrix<N> operator-(const Matrix<N> a, const Matrix<N>& b) {
- return subtract(a, b);
-}
+ template <size_t N>
+ inline Matrix<N> operator-(const Matrix<N> a, const Matrix<N>& b) {
+ return subtract(a, b);
+ }
-template <size_t N>
-inline Matrix<N> multiply(const double a, const Matrix<N>& m) {
- Matrix<N> rtn;
- rtn._matrix = a * m._matrix;
- return rtn;
-}
+ template <size_t N>
+ inline Matrix<N> multiply(const double a, const Matrix<N>& m) {
+ Matrix<N> rtn;
+ rtn._matrix = a * m._matrix;
+ return rtn;
+ }
-template <size_t N>
-inline Matrix<N> multiply(const Matrix<N>& m, const double a) {
- return multiply(a, m);
-}
+ template <size_t N>
+ inline Matrix<N> multiply(const Matrix<N>& m, const double a) {
+ return multiply(a, m);
+ }
-template <size_t N>
-inline Matrix<N> divide(const Matrix<N>& m, const double a) {
- return multiply(1/a, m);
-}
+ template <size_t N>
+ inline Matrix<N> divide(const Matrix<N>& m, const double a) {
+ return multiply(1/a, m);
+ }
-template <size_t N>
-inline Matrix<N> operator*(const double a, const Matrix<N>& m) {
- return multiply(a, m);
-}
+ template <size_t N>
+ inline Matrix<N> operator*(const double a, const Matrix<N>& m) {
+ return multiply(a, m);
+ }
-template <size_t N>
-inline Matrix<N> operator*(const Matrix<N>& m, const double a) {
- return multiply(a, m);
-}
+ template <size_t N>
+ inline Matrix<N> operator*(const Matrix<N>& m, const double a) {
+ return multiply(a, m);
+ }
-template <size_t N>
-inline Matrix<N> multiply(const Matrix<N>& a, const Matrix<N>& b) {
- Matrix<N> tmp;
- tmp._matrix = a._matrix * b._matrix;
- return tmp;
-}
+ template <size_t N>
+ inline Matrix<N> multiply(const Matrix<N>& a, const Matrix<N>& b) {
+ Matrix<N> tmp;
+ tmp._matrix = a._matrix * b._matrix;
+ return tmp;
+ }
-template <size_t N>
-inline Matrix<N> operator*(const Matrix<N>& a, const Matrix<N>& b) {
- return multiply(a, b);
-}
+ template <size_t N>
+ inline Matrix<N> operator*(const Matrix<N>& a, const Matrix<N>& b) {
+ return multiply(a, b);
+ }
-template <size_t N>
-inline Vector<N> multiply(const Matrix<N>& a, const Vector<N>& b) {
- Vector<N> tmp;
- tmp._vec = a._matrix * b._vec;
- return tmp;
-}
+ template <size_t N>
+ inline Vector<N> multiply(const Matrix<N>& a, const Vector<N>& b) {
+ Vector<N> tmp;
+ tmp._vec = a._matrix * b._vec;
+ return tmp;
+ }
-template <size_t N>
-inline Vector<N> operator*(const Matrix<N>& a, const Vector<N>& b) {
- return multiply(a, b);
-}
+ template <size_t N>
+ inline Vector<N> operator*(const Matrix<N>& a, const Vector<N>& b) {
+ return multiply(a, b);
+ }
-template <size_t N>
-inline Matrix<N> transpose(const Matrix<N>& m) {
- // Matrix<N> tmp;
- // for (size_t i = 0; i < N; ++i) {
- // for (size_t j = 0; j < N; ++j) {
- // tmp.set(i, j, m.get(j, i));
- // }
- // }
- // return tmp;
- return m.transpose();
-}
+ template <size_t N>
+ inline Matrix<N> transpose(const Matrix<N>& m) {
+ // Matrix<N> tmp;
+ // for (size_t i = 0; i < N; ++i) {
+ // for (size_t j = 0; j < N; ++j) {
+ // tmp.set(i, j, m.get(j, i));
+ // }
+ // }
+ // return tmp;
+ return m.transpose();
+ }
-template <size_t N>
-inline Matrix<N> inverse(const Matrix<N>& m) {
- return m.inverse();
-}
+ template <size_t N>
+ inline Matrix<N> inverse(const Matrix<N>& m) {
+ return m.inverse();
+ }
-template <size_t N>
-inline double det(const Matrix<N>& m) {
- return m.determinant();
-}
+ template <size_t N>
+ inline double det(const Matrix<N>& m) {
+ return m.determinant();
+ }
-template <size_t N>
-inline double trace(const Matrix<N>& m) {
- return m.trace();
-}
+ template <size_t N>
+ inline double trace(const Matrix<N>& m) {
+ return m.trace();
+ }
-/////////////////////////////////
+ /////////////////////////////////
-/// Make string representation
-template <size_t N>
-inline string toString(const Matrix<N>& m) {
- ostringstream ss;
- ss << "[ ";
- for (size_t i = 0; i < m.size(); ++i) {
- ss << "( ";
- for (size_t j = 0; j < m.size(); ++j) {
- const double e = m.get(i, j);
- ss << (Rivet::isZero(e) ? 0.0 : e) << " ";
+ /// Make string representation
+ template <size_t N>
+ inline string toString(const Matrix<N>& m) {
+ ostringstream ss;
+ ss << "[ ";
+ for (size_t i = 0; i < m.size(); ++i) {
+ ss << "( ";
+ for (size_t j = 0; j < m.size(); ++j) {
+ const double e = m.get(i, j);
+ ss << (Rivet::isZero(e) ? 0.0 : e) << " ";
+ }
+ ss << ") ";
}
- ss << ") ";
+ ss << "]";
+ return ss.str();
}
- ss << "]";
- return ss.str();
-}
-/// Stream out string representation
-template <size_t N>
-inline ostream& operator<<(std::ostream& out, const Matrix<N>& m) {
- out << toString(m);
- return out;
-}
+ /// Stream out string representation
+ template <size_t N>
+ inline ostream& operator<<(std::ostream& out, const Matrix<N>& m) {
+ out << toString(m);
+ return out;
+ }
-/////////////////////////////////////////////////
+ /////////////////////////////////////////////////
-/// Compare two matrices by index, allowing for numerical precision
-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;
+ /// Compare two matrices by index, allowing for numerical precision
+ 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;
}
- return true;
-}
-/// External form of numerically safe nullness check
-template <size_t N>
-inline bool isZero(const Matrix<N>& m, double tolerance=1E-5) {
- return m.isZero(tolerance);
-}
+ /// External form of numerically safe nullness check
+ template <size_t N>
+ inline bool isZero(const Matrix<N>& m, double tolerance=1E-5) {
+ return m.isZero(tolerance);
+ }
}
-
#endif
Modified: trunk/include/Rivet/Math/StdHeader.hh
==============================================================================
--- trunk/include/Rivet/Math/StdHeader.hh Wed Nov 4 00:12:46 2009 (r1991)
+++ trunk/include/Rivet/Math/StdHeader.hh Wed Nov 4 00:46:05 2009 (r1992)
@@ -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 Wed Nov 4 00:12:46 2009 (r1991)
+++ trunk/include/Rivet/Math/VectorN.hh Wed Nov 4 00:46:05 2009 (r1992)
@@ -4,7 +4,7 @@
#include "Rivet/Math/MathHeader.hh"
#include "Rivet/Math/MathUtils.hh"
-#include <Eigen/Core>
+#include "Rivet/Math/eigen/vector.h"
namespace Rivet {
@@ -25,22 +25,21 @@
friend Vector<M> multiply(const Matrix<M>& a, const Vector<M>& b);
public:
- Vector() { _vec.setZero(); }
+ Vector() { _vec.loadZero(); }
Vector(const Vector<N>& other)
: _vec(other._vec) { }
- /// Get vector elements
- const double get(const size_t index) const {
+ 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);
}
@@ -128,13 +127,12 @@
if (index >= N) {
throw std::runtime_error("Tried to access an invalid vector index.");
} else {
- return _vec[index];
+ return _vec(index);
}
}
/// Vector
- typedef Eigen::Matrix<double,N,1> EVector;
- EVector _vec;
+ Eigen::Vector<double,N> _vec;
};
Modified: trunk/pyext/setup.py.in
==============================================================================
--- trunk/pyext/setup.py.in Wed Nov 4 00:12:46 2009 (r1991)
+++ trunk/pyext/setup.py.in Wed Nov 4 00:46:05 2009 (r1992)
@@ -15,7 +15,8 @@
ext = Extension('_rivet',
['@srcdir@/rivet_wrap.cc'],
define_macros = [("SWIG_TYPE_TABLE", "hepmccompat")],
- include_dirs=[incdir, os.path.join(incdir, 'eigen2'), '@HEPMCINCPATH@', '@BOOSTINCPATH@', '@GSLINCPATH@'],
+ #include_dirs=[incdir, os.path.join(incdir, 'eigen2'), '@HEPMCINCPATH@', '@BOOSTINCPATH@', '@GSLINCPATH@'],
+ include_dirs=[incdir, '@HEPMCINCPATH@', '@BOOSTINCPATH@', '@GSLINCPATH@'],
#should replace '.libs' -> os.path.join(srcdir,'@LT_OBJDIR@'), but doesn't work
library_dirs=[srcdir, os.path.join(srcdir,'.libs'), '@HEPMCLIBPATH@', '@GSLLIBPATH@'],
libraries=['HepMC', 'Rivet'])
More information about the Rivet-svn
mailing list
|