00001 #ifndef RIVET_MATH_VECTORN
00002 #define RIVET_MATH_VECTORN
00003
00004 #include "Rivet/Math/MathHeader.hh"
00005 #include "Rivet/Math/MathUtils.hh"
00006
00007 #include "Rivet/Math/eigen/vector.h"
00008
00009 namespace Rivet {
00010
00011
00012 template <size_t N>
00013 class Vector;
00014 template <size_t N>
00015 class Matrix;
00016
00017 template <size_t N>
00018 Vector<N> multiply(const Matrix<N>& a, const Vector<N>& b);
00019
00020
00021
00022 template <size_t N>
00023 class Vector {
00024 template <size_t M>
00025 friend Vector<M> multiply(const Matrix<M>& a, const Vector<M>& b);
00026
00027 public:
00028 Vector() { _vec.loadZero(); }
00029
00030 Vector(const Vector<N>& other)
00031 : _vec(other._vec) { }
00032
00033 const double& get(const size_t index) const {
00034 if (index >= N) {
00035 throw std::runtime_error("Tried to access an invalid vector index.");
00036 } else {
00037 return _vec(index);
00038 }
00039 }
00040
00041
00042 const double& operator[](const size_t index) const {
00043 return get(index);
00044 }
00045
00046
00047 double& operator[](const size_t index) {
00048 return get(index);
00049 }
00050
00051
00052 Vector<N>& set(const size_t index, const double value) {
00053 if (index >= N) {
00054 throw std::runtime_error("Tried to access an invalid vector index.");
00055 } else {
00056 _vec[index] = value;
00057 }
00058 return *this;
00059 }
00060
00061
00062 size_t size() const {
00063 return N;
00064 }
00065
00066
00067 bool isZero(double tolerance=1E-5) const {
00068 for (size_t i=0; i < N; ++i) {
00069 if (! Rivet::isZero(_vec[i], tolerance) ) return false;
00070 }
00071 return true;
00072 }
00073
00074
00075
00076 double mod2() const {
00077 double mod2 = 0.0;
00078 for (size_t i = 0; i < size(); ++i) {
00079 const double element = get(i);
00080 mod2 += element*element;
00081 }
00082 return mod2;
00083 }
00084
00085
00086
00087 double mod() const {
00088 const double norm = mod2();
00089 assert(norm >= 0);
00090 return sqrt(norm);
00091 }
00092
00093
00094 Vector<N> operator-() const {
00095 Vector<N> rtn;
00096 rtn._vec = -_vec;
00097 return rtn;
00098 }
00099
00100 bool operator==(const Vector<N>& a) const {
00101 return _vec == a._vec;
00102 }
00103
00104 bool operator!=(const Vector<N>& a) const {
00105 return _vec != a._vec;
00106 }
00107
00108 bool operator<(const Vector<N>& a) const {
00109 return _vec < a._vec;
00110 }
00111
00112 bool operator<=(const Vector<N>& a) const {
00113 return _vec <= a._vec;
00114 }
00115
00116 bool operator>(const Vector<N>& a) const {
00117 return _vec > a._vec;
00118 }
00119
00120 bool operator>=(const Vector<N>& a) const {
00121 return _vec >= a._vec;
00122 }
00123
00124
00125 protected:
00126 double& get(const size_t index) {
00127 if (index >= N) {
00128 throw std::runtime_error("Tried to access an invalid vector index.");
00129 } else {
00130 return _vec(index);
00131 }
00132 }
00133
00134
00135 Eigen::Vector<double,N> _vec;
00136 };
00137
00138
00139
00140
00141
00142 template <size_t N>
00143 inline double mod2(const Vector<N>& v) {
00144 return v.mod2();
00145 }
00146
00147
00148
00149 template <size_t N>
00150 inline double mod(const Vector<N>& v) {
00151 return v.mod();
00152 }
00153
00154
00155
00156
00157
00158
00159 template <size_t N>
00160 inline const string toString(const Vector<N>& v) {
00161 ostringstream out;
00162 out << "(";
00163 for (size_t i = 0; i < v.size(); ++i) {
00164 out << (fabs(v[i]) < 1E-30 ? 0.0 : v[i]);
00165 if (i < v.size()-1) out << ", ";
00166 }
00167 out << ")";
00168 return out.str();
00169 }
00170
00171
00172 template <size_t N>
00173 inline std::ostream& operator<<(std::ostream& out, const Vector<N>& v) {
00174 out << toString(v);
00175 return out;
00176 }
00177
00178
00179
00180
00181
00182
00183 template <size_t N>
00184 inline bool fuzzyEquals(const Vector<N>& va, const Vector<N>& vb, double tolerance=1E-5) {
00185 for (size_t i = 0; i < N; ++i) {
00186 const double a = va.get(i);
00187 const double b = vb.get(i);
00188 if (!Rivet::fuzzyEquals(a, b, tolerance)) return false;
00189 }
00190 return true;
00191 }
00192
00193
00194
00195 template <size_t N>
00196 inline bool isZero(const Vector<N>& v, double tolerance=1E-5) {
00197 return v.isZero(tolerance);
00198 }
00199
00200
00201 }
00202
00203 #endif