Rivet/Math/VectorN.hh

0001 #ifndef RIVET_MATH_VECTORN
0002 #define RIVET_MATH_VECTORN
0003
0004 #include "Rivet/Math/MathConstants.hh"
0005 #include "Rivet/Math/MathUtils.hh"
0006
0007 #include "Rivet/Math/eigen3/Dense"
0008
0009 namespace Rivet {
0010
0011
0012   template <size_t N>
0013   class Vector;
0014   template <size_t N>
0015   class Matrix;
0016
0017   template <size_t N>
0018   Vector<N> multiply(const Matrix<N>& a, const Vector<N>& b);
0019
0020
0021   /// A minimal base class for \f$ N \f$-dimensional vectors.
0022   template <size_t N>
0023   class Vector {
0024
0025     template <size_t M>
0026     friend Vector<M> multiply(const Matrix<M>& a, const Vector<M>& b);
0027
0028
0029   public:
0030
0031     Vector() : _vec(EVector::Zero()) { }
0032
0033     const double& get(const size_t index) const {
0034       if (index >= N) {
0035         throw std::runtime_error("Tried to access an invalid vector index.");
0036       } else {
0037         return _vec(index);
0038       }
0039     }
0040
0041     double& get(const size_t index) {
0042       if (index >= N) {
0043         throw std::runtime_error("Tried to access an invalid vector index.");
0044       } else {
0045         return _vec(index);
0046       }
0047     }
0048
0049     /// Direct access to vector elements by index.
0050     const double& operator[](const size_t index) const {
0051       return get(index);
0052     }
0053
0054     /// Direct access to vector elements by index.
0055     double& operator[](const size_t index) {
0056       return get(index);
0057     }
0058
0059     /// Set indexed value
0060     Vector<N>& set(const size_t index, const double value) {
0061       if (index >= N) {
0062         throw std::runtime_error("Tried to access an invalid vector index.");
0063       } else {
0064         _vec[index] = value;
0065       }
0066       return *this;
0067     }
0068
0069     /// Vector dimensionality
0070     constexpr size_t size() const {
0071       return N;
0072     }
0073
0074     /// Check for nullness, allowing for numerical precision
0075     bool isZero(double tolerance=1E-5) const {
0076       for (size_t i=0; i < N; ++i) {
0077         if (! Rivet::isZero(_vec[i], tolerance) ) return false;
0078       }
0079       return true;
0080     }
0081
0082     /// @brief Calculate the modulus-squared of a vector.
0083     /// \f$ \sum_{i=1}^N x_i^2 \f$.
0084     double mod2() const {
0085       double mod2 = 0.0;
0086       for (size_t i = 0; i < size(); ++i) {
0087         const double element = get(i);
0088         mod2 += element*element;
0089       }
0090       return mod2;
0091     }
0092
0093     /// @brief Calculate the modulus of a vector.
0094     /// \f$ \sqrt{\sum_{i=1}^N x_i^2} \f$.
0095     double mod() const {
0096       const double norm = mod2();
0097       //assert(norm >= 0); //< *should* be impossible
0098       return sqrt(norm);
0099     }
0100
0101     /// Invert the vector
0102     Vector<N> operator-() const {
0103       Vector<N> rtn;
0104       rtn._vec = -_vec;
0105       return rtn;
0106     }
0107
0108     bool operator==(const Vector<N>& a) const {
0109       return _vec == a._vec;
0110     }
0111
0112     bool operator!=(const Vector<N>& a) const {
0113       return _vec != a._vec;
0114     }
0115
0116     // bool operator<(const Vector<N>& a) const {
0117     //   return _vec < a._vec;
0118     // }
0119
0120     // bool operator<=(const Vector<N>& a) const {
0121     //   return _vec <= a._vec;
0122     // }
0123
0124     // bool operator>(const Vector<N>& a) const {
0125     //   return _vec > a._vec;
0126     // }
0127
0128     // bool operator>=(const Vector<N>& a) const {
0129     //   return _vec >= a._vec;
0130     // }
0131
0132     /// Vector
0133     using EVector = RivetEigen::Matrix<double,N,1>;
0134     EVector _vec;
0135
0136   };
0137
0138
0139   /////////////////////////////////////////////////
0140
0141
0142   /// @name String representations of vectors
0143   /// @{
0144
0145   /// Make string representation
0146   template <size_t N>
0147   inline const string toString(const Vector<N>& v) {
0148     std::ostringstream out;
0149     out << "(";
0150     for (size_t i = 0; i < v.size(); ++i) {
0151       out << (fabs(v[i]) < 1E-30 ? 0.0 : v[i]);
0152       if (i < v.size()-1) out << ", ";
0153     }
0154     out << ")";
0155     return out.str();
0156   }
0157
0158   /// Stream out string representation
0159   template <size_t N>
0160   inline std::ostream& operator<<(std::ostream& out, const Vector<N>& v) {
0161     out << toString(v);
0162     return out;
0163   }
0164
0165   /// @}
0166
0167
0168   /////////////////////////////////////////////////
0169
0170
0171   /// Compare two vectors by index, allowing for numerical precision
0172   template <size_t N>
0173   inline bool fuzzyEquals(const Vector<N>& va, const Vector<N>& vb, double tolerance=1E-5) {
0174     for (size_t i = 0; i < N; ++i) {
0175       const double a = va.get(i);
0176       const double b = vb.get(i);
0177       if (!Rivet::fuzzyEquals(a, b, tolerance)) return false;
0178     }
0179     return true;
0180   }
0181
0182
0183   /// External form of numerically safe nullness check
0184   template <size_t N>
0185   inline bool isZero(const Vector<N>& v, double tolerance=1E-5) {
0186     return v.isZero(tolerance);
0187   }
0188
0189
0190 }
0191
0192 #endif